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Vorwort 


Kaum ein technisches Gerät hat sich so explosionsartig entwickelt 
wie der elektronische Taschenrechner. Innerhalb von fünf Jahren 
wuchs er von einem Investitionsgut für einen kleinen Kreis von 
Spezialisten zu einem Massenprodukt heran, das, in Kaufhäusern 
angeboten, einen völlig neuen und bisher unbekennten Narkt be- 
dient. Doch ihre volle Ausnutzung und Bedeutung erlangen diese 
Geräte erst durch die Kenntnis der zugehörigen numerischen Ver- 
fahren. Nur langsam folgt hier die notwendige Literatur der all- 
zuraschen Entwicklung. So versucht dieses Buch, eine Lücke kurz- 
fristig zu schließen. Dabei stellte sich die reizvolle Aufgabe, 
eine Brücke über mehrere Jahrhunderte hinweg zu schlagen, zwi- 
schen alten, sehr einfachen Rechensystemen und der neuen, elek- 
tronischen Rechentechnik, 

Das Buch richtet sich an Schüler der gymnasialen Oberstufe, an 
Lehrende, Lernende und Praktizierende der Mathematik, Physik und 
Ingenieurwissenschaften und möchte bei der Lösung ihrer Probleme 
helfen. Deshalb wechselt bei den Zahlenbeispielen Bekanntes mit 
weniger Bekanntem und Neuem ab. Die Beispiele wurden auf einem 
Taschenrechner HP 67 der Firma Hewlett Packard gerechnet. Die 
Konzeption dieses Rechners eignet sich besonders zur Lösung par- 
tieller Differentialgleichungen. 

Um den Text durch theoretische Erörterungen nicht zu belasten, 
wurden die wichtigsten mathematischen Existenzbedingungen - 

2.B. Lipschitzbedingung - stillschweigend vorausgesetzt. 

Dem Verlag danke ich für seine freundliche Hilfe bei der Arbeit 
an diesem Buch. 


Bremen, im Frühjahr 1978 Genhayd Vene 


Einleitung 


von Naturvorgängen die wohl am häufigsten anzutreffende Glei- 
chungsform, weil sich die physikalischen Gesetze am deutlichstenr 
in ihnen ausdrücken lassen. So wird z.B. das Einschalten einer 
Spule mit der Induktivitäöt L über einen Widerstand R durch die 
Differentialgleichung 


ee si 
u=iR+L at 


beschrieben. Um so erfolgreicher sträuben sich aber Differential- 
gleichungen gegen ihre formelmäßige Auflösung. Ja, das Auffinden 
einer analytischen Lösung mit Hilfe bekannter Funktionen, wie x”, 
In x oder sin x „ gehört zu den seltenen Ausnahmen. Hier helfen 
nur zeichnerische oder numerische Näherungslösungen weiter, wobei 
den numerischen Verfahren der Vorzug zu geben ist, wegen ihrer 
Vielseitigkeit und weil die Lösung der Differentialgleichung mit 
gewünschter Genauigkeit angegeben werden kann. 

Die numerischen Verfahren beinhalten aber einen großen Rechenauf- 
wand. Erst durch die Erstellung von Rechenautomaten, mit der Mög- 
lichkeit, sehr viele Rechenoperationen in kürzester Zeit auszufüh- 
ren, wurden die numerischen Verfahren für den Praktiker inter- 
essant. Mit der Weiterentwicklung zu programmierbaren Rechnern 
im Taschenformat eröffnet sich nunmehr dem Interessierten die 
Möglichkeit, diesen reizvollen Zweig der Naturbeschreibung durch 


R 


Abb. 1 Einschalten einer Spule 
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die Lösung vorliegender Differentialgleichungen auch außerhalb 
von Rechenzentren kennenzulernen. 

Man unterscheidet gewöhnliche und partielle Differentialgleichun- 
gen. Bei den gewöhnlichen Differentialgleichungen kommen nur Ab- 
leitungen (Differentialquotienten) nach einer Veränderlichen vor, 
z.B. der Zeit t oder einer Auslenkung x. Bei partiellen Differen- 
tialgleichungen finden sich dagegen Ableitungen nach mehreren Ver- 
änderlichen. Für beide Arten werden numerische Verfahren beschrie- 
ben, die sich für programmierbare Taschenrechner mit ihrer be- 
schränkten Anzahl von Speicherplätzen besonders gut eignen. Die 
untere Grenze des notwendigen Programmierumfanges liegt bei etwa 
250 Programmschritten (ROM Read Only Memory) und 25 Datenspei- 
chern (RAM Random Access Memory); außerdem müssen indirektes 
Adressieren, Programmverzweigungen (Decision) nach Vergleichs- 
operationen und Unterprogramme (Subroutine) möglich sein. 


Teil A 


Gewöhnliche Differentialgleichungen 


1.  Differentialgleichungen mit Anfangswerten 


11  Differentialgleichungen erster Ordnung 


Eine Differentialgleichung erster Ordnung ist eine Gleichung zwi- 
schen den Veränderlichen x und y und der ersten Ableitung y' nach 
x. Sie sei umstellbar in die explizite Form 


y'=f(x,y) - 


Gesucht wird die Kurve y(x), die die Differentialgleichung er- 
füllt und durch einen vorgegebenen Anfangspunkt x» Yo geht. 
Diese Kurve soll punktweise, jeweils im Abstand einer in x-Rich- 
tung gewählten Schrittweite h errechnet werden (Abb.2). 

Aus der Vielzahl der bekannten numerischen Lösungswege eignet 
sich das Verfahren nach RUNGE und KUTTA für programmierbare Ta- 
schenrechner besonders gut, aus folgenden Gründen: 


Es ist zuverlässig, arbeitet rein mechanisch und ist einfach zu 
programmieren, 

Es arbeitet ohne Iteration und Anlaufrechnung. 

Es zeichnet sich bei richtiger Schrittwahl h durch eine hohe Ge- 
nauigkeit aus. 

Es erlaubt im Bedarfsfall das Abändern der Schrittweite h nach 
jedem berechneten Punkt der Kurve y(x). 


[e} 


Abb. 2 Punktweise Berechnung einer Kurve 
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Es gestattet dadurch, bei geringem Mehraufwand, eine automati- 
sche Steuerung zur richtigen Schrittweite hin. 

Es läßt ohne weiteres die Erweiterung auf Differentialgleichun- 
gen zweiter und höherer Ordnung zu. 


Auf die mathematische Begründung dieses Verfahrens wird hier ver- 
zichtet, weil sie einen erheblichen Aufwand erfordert. Man kann 
sie aber nachlesen, z.B. in dem Buch von C. Runge und H. König, 
Vorlesungen über numerisches Rechnen, Berlin 1924. 


Beschreibung des Verfahrens 


Das Verfahren nach RUNGE-KUTTA beruht auf einem dreimaligen Sich- 
Vortasten längs verschiedener Parabelbögen, beginnend mit dem An- 
fangspunkt X Yo» und zwar zweimal un bis zur Schrittmitte 
und einmal um h bis zum Schrittende (Abb.3). Hierbei werden an 
den Schrittstellen S] vier verschiedene Zuwächse k] (I= 1 bis 4) 
in y-Richtung berechnet, aus denen ein Mittelwert k gebildet wird. 
Mit h und k findet man das Wertepaar des nächsten Punktes 


x =Xx, + h yı =, + k. 


Dieses dient als neuer Anfangspunkt für die nächste Schrittbe- 
rechnung usw. 


Abb. 3 RUNGE-KUTTA-Verfahren 
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Das Gleichungssystem für dieses Vorgehen wurde für die Zahlen- 
rechnung auf eine besonders günstige Form gebracht und schreibt 


sich an einer beliebigen Stelle X r In? 
Zuwachs 
K- =2.f(s,) mit S1: X > In 
k, -$-f(s,) 52: %, +; „+ KK (1.1) 
a, . a 
K, =57 f(s,) SE nt mr k, 
k A.rls ) s,: x +h y_+2k 
472 4 4° n . n 3 >? 
Mittelwert 
4 
k = z.( kı + 2:k, + 2:k, 4 ky Je; (1.2) 


Wertepaar des nächsten Punktes 


Kart = X +h Yn+1 = Ya + K . (1.3) 


Programmablaufplan 


Wie wir gesehen haben, besteht das RUNGE-KUTTA-Verfahren aus ei- 
ner Reihe von Rechenoperationen und logischen Entscheidungen, die 
das Wertepaar Kn9 Ya eines Punktes der Lösungskurve mit dem Wer- 
tepaar Karl r Int des nächsten Punktes verknüpft. Ein wertvol- 
les Hilfsmittel für die Umsetzung dieser einzelnen Rechenschritte 
in das Programm eines Taschenrechners ist der Programmablaufplan 
oder das Flußdiagramm. In ihm wird der Lösungsablauf graphisch 


dargestellt. So zeigt der Programmablaufplan insbesondere die 


Ein- und Ausgabe von Zahlen, 
Rechenoperationen, 
logischen Entscheidungen, 


und zwar in der Reihenfolge und Verknüpfung, in der sie bei der 
Zahlenrechnung ablaufen müssen. An Hand eines Flußdiagramnes, 
läßt sich das Programm für einen Taschenrechner Schritt für 
Schritt erstellen. 

Für die graphische Darstellung des Lösungsablaufes wurden im 
DIN-Blatt 66001 verschiedene Symbole vereinbart, von denen neun 
in diesem Buch verwendet werden: 
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Jo-+! OILPL 


Das Rechteck hebt die Stellen hervor, an denen 
Rechenoperationen ausgeführt werden. 


Ein Rhombus kennzeichnet eine Verzweigungsstelle. 


Dieses Zeichen bedeutet das Abarbeiten eines 
Unterprogrammes, ehe im Hauptprogramm weiter- 
gegangen wird. 


Stellen, an denen Daten ein- oder ausgegeben 
werden, kennzeichnet das Parallelogramn. 


Der Strich bedeutet eine Ablauflinie. Vorzugs- 
richtung ist von oben nach unten oder von links 
nach rechts. 


So wird eine Zusammenführung gezeichnet. Zwei 
sich kreuzende Ablauflinien bedeuten keine Zu- 
sammenführung. 


Der Kreis steht für eine Übergangsstelle. Zu- 
sammengehörige Übergangsstellen müssen gleiche 
Bezeichnungen haben. 


Mit diesem Zeichen werden Anfang und Ende, auch 
Zwischenhalte einer vollständigen Rechenarbeit 
markiert. 


Dieses Zeichen wird für Bemerkungen an irgend- 
ein Symbol angefügt. 
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Erhöhe 


I um 1 


Abb. A Differentialgleichung erster Ordnung 
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Den Programmablaufplan für eine Differentialgleichung erster Ord- 
nung mit Anfangswerten zeigt die Abb.4. Um mit dem gleichen Pro- 
gramm die unterschiedlichsten Differentialgleichungen 


y' = f(x,y) 


lösen zu können, ist im Ablaufplan für die Funktion f(x,y) ein 
Unterprogramm vorgesehen. Ist das Hauptprogramm einmal erstellt, 
so genügt es, zur Lösung einer vorliegenden Differentialgleichung 
die Funktion f(x,y) als Unterprogramm zu schreiben und das Haupt- 
programm mit den Anfangswerten x Y, Zu starten. 


1.1.1 Beispiel. Oberfläche von Flüssigkeitsströmung 


Die Oberfläche einer gestauten Wasserströmung im längsschnitt 
(Abb.5) wird durch die Differentialgleichung 


1b -8)] 


beschrieben, sofern ein seichtes Bett von parabolischem Quer- 
schnitt vorliegt. Gesucht wird die Wasserhöhe y im Abstand von 
der Staumauer x für die folgenden Zahlenwerte: 


Y = 0,0873 kleines Gefälle im Bogenmaß 

= 1 m Wasserhöhe in größerer Entfernung x 

1,5m für x = O0 m Anfangswerte 

= 1 m Schrittweite für das RUNGE-KUTTA-Verfahren. 


H 
Yo 
h 


Ergebnis: 


OVOVO@-HINUNPBWND-0O 53 4% 
pe 
D 
\n 
KW 


_ 


1.1 Differentialgleichungen erster Ordnung 19 


4 8 12 m 


Abb. 5 Flüssigkeitsströmung 


Die Zahlen wurden auf einem HP 67 mit dem folgenden Programm ge- 


rechnet: 


Hauptprogramm Unterprogramm f(x,y) 
LBL A GSB O0 RCL 3 GSB B LBL B 
RCL 1 2 STO+2 RCL O RCL 6 
STO A x GTO A x RCLB 
—X- RCLO STO+3 s 
RCL 2 STO+1 LBL 0 STO+3 4 
STO B R4 RCL 2 RTN x 
=I= GSB 0 + Y 
GSB B STO-3 STO a e 
RCLO 3 RCL 
x 570: 3 RCL 0 en 
STo 3 RCL © + ih 
GSB 0 STO+1 STO A 

y 


Yan Grobrechnung 


Yn+h Feinrechmung 


Abb. 6 Fehlerabschätzung 
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Speicherplan 


A 
x 


Wertepaar 


belegt frei 


o.gla 
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ergleich 
mit € 


B 


nach (1.6) 


Verdoppeln 


Rückruf 
Healbieren 


8 
A 
25 
>» 
a 
Ag 
vo © 
m» 


Abb. 7 Differentialgleichung erster Ordnung mit Korrektur und 
Schrittsteuerung 


1.1.2 Schrittsteuerung und Korrektur 


Das RUNGE-KUTTA-Verfahren ist hinsichtlich seiner Genauigkeit ein 
Verfahren vierter Ordnung, d.h. der Fehler ist von der Ordnung n°. 
Bei hinreichend kleiner Schrittweite wird daher der Benutzer 
durch eine hohe Genauigkeit belohnt. Umgekehrt verhält es sich 
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aber bei zu großer Schrittweite, Die richtige Schrittbemessung 
ist bei diesem Verfahren von entscheidender Bedeutung. Hier liegt 
eine Quelle der Unsicherheit, zumal sich eine Fehlerrechnung nur 
mit großer Schwierigkeit ausführen läßt. 

Dagegen gibt es für ein angenähertes Abschätzen des Fehlers fol- 
genden einfachen Weg: 

Die Rechnung wird mit doppelter Schrittweite 2h wiederholt und 
das Ergebnis Y>h dieser groben Rechnung mit der zuvor durchge- 
führten Feinrechnung Yn+h verglichen. Dann läßt sich angenähert 
sagen, daß der Fehler bei der Schrittverdoppelung auf den 2?- 32- 
fachen Wert angestiegen ist. Da gleichzeitig zum Durchlaufen von 
2h zwei Schritte h + h der Feinrechnung nötig sind, wird sich 
hier der Fehler nahezu verdoppelt haben. So kann man den Fehler 
der Grobrechnung auf das 16fache der Feinrechnung schätzen. Die 
Differenz zwischen Grob- und Feinrechnung wird somit etwa das 
15fache des Fehlers in der Feinrechnung betragen (Abb.6). 

Diese Fehlerabschätzung läßt sich sowohl zur automatischen 
Schrittsteuerung als auch zur Korrektur der Feinrechnung heran- 
ziehen. 


Korrektur der Feinrechnung 
1 
8 = 15.( Yn+h - Yan)" (1.4) 
Korrigierter Wert 


(1.5) 


1 
ll 


x +h+h. 


für n 


{ai 
u 


Nach der Wahl einer Genauigkeitsschranke e, z.B. vier Stellen 
nach dem Komma, hat sich für die Schrittsteuerung folgende Fest- 
legung als brauchbar herausgestellt /9/ : 


IE < |ö] (a) beibehalten von h 
2 >18 (b) verdoppeln von h (1.6) 
10€e< |ö| (ec) halbieren vonh. 


Im Programmablaufplan (Abb.7) für eine Rechnung mit Schrittsteu- 
erung und Korrektur ist zur Vereinfachung die Berechnung des Mit- 
telwertes k als Unterprogramm gezeichnet. 
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1,1.3 Beispiel. Einschalten einer Spule mit Eisenkern 


Bei einer Spule mit Eisenkern hängt die Induktivität L vom Sätti- 
gungsgrad des magnetischen Flusses $ ab, der vom Strom i bestimmt 
wird. Dann erhält die Differentialgleichung zur Abb. 1 die Form 


u=-iR+ moi) . N Windungszahl der Spule 
Zur analytischen Nachbildung der Magnetisierungskurve P in einen 
magnetischen Kreis aus Eisen und luftspalt ist die Funktion 


F(z) = coth(z) -4 
gut geeignet, deren Kurvenverlauf die Abb. 8 zeigt. Für die wei- 


tere Rechnung wird auch noch ihre Ableitung benötigt 


L: VoBE0per DERBER L 
4 „2 sinh’(z) 
gar 


F(z) = eoth(z) - 
0.5 


Abb. 8 Nachbildung der Eisen-Magnetisierungskurve 
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Mit $ = Po F(z) und z= @-i erhält man eine veränderliche In- 
duktivität 


L(i) = N®,« 5 
mit der sich die Differentialgleichung schreibt 
Be ).& 
u=i'R+ L(i) at 


Umgestellt für das RUNGE-KUTTA-Verfahren 


Die Induktivität L(i) beinhaltet selber eine Funktion, die nicht 
ganz einfach aufgebaut ist. Für die Rechnung nach RUNGE-KUTTA 
bedeuten nichtkonstante Faktoren aber keine zusätzliche Erschwer- 
nis. 


Die Lösung wird für folgende Zahlenwerte gesucht: 


NG,% 1 1 
E12; ——: 0,3326=55; a=0 g und a=5%; 
i, =0A für t, = 0 s Anfangswerte; h = 0,5 s erste Schrittweite. 
Ergebnis: 


2 
u 
oO 
R 
) 
\n 


0,00 
0,19 
0,65 
0,97 
0,99 
1,00 


NIJANW-O U cd 


Aus den Zahlenreihen für die Zeit + ist die automatische Schritt- 
steuerung deutlich ablesbar. Aus der Abb. 9 ist ersichtlich, daß 

der Einschaltstrom i bei höherem Sättigungsgrad schneller seinen 

Endwert erreicht. 

Die Zahlen wurden auf einem HP 67 mit dem folgenden Programm ge- 

rechnet: 


1.1 Differentialgleichungen ersıer Ordrurg 


Abb. 9 Einschaltstrom einer Spule mit Eisenkern 


Hauptprogramm 
IBLA LBL 2 
3 RCL 1 
ST I STO A 
RCL 1 RCL 2 
STO A STO B 
—- GSB B 
RCL 2 RCLO 
STO B x 

STO 4 STO 3 
—- GSB O 
GTO 2 GSB O0 
IBL 1 2 

RCL O x 

2 RCLO 
x STO+1 
STO-1 R} 
STO-1 GSB O0 
STO O0 STO-3 
RCL 4 3 
RCL 2 STO: 3 
STO 4 RCL Oo 
Rh STO+1 
STO 2 RCL 3 


STO+r2 


x 


x=0 

GTO b 

ENTER 
1/x 
2 


x 


aY 


e 


ENTER 
1/x 


2 


x 


x 


1/x 


RCL 
x 
RCL 


Unterprogramn f(x,y) 


LBL B 
RAD 

RCL B 
RCL5 


7 
6 


8 
B 


A..Wo0c 
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Speicherplan 
A B C D E 

Wertepear x y frei frei frei 

0 2 3 4 5 6 7 3 9 

h u ; 
> x, AR belegt belest % 0,3326 5 s R frei 
10 bis 19 

frei 


1.2 Gekoppelte Differentialgleichungen erster Ordnung, 
Differentialgleichungen zweiter Ordnung 


1.2.1 Problemstellung und Lösungsweg 


Häufig wird ein physikalischer Vorgang nicht durch eine einzige 
Differentialgleichung beschrieben, sondern durch ein System von 
zwei oder mehreren Differentialgleichungen, die miteinander zu- 
sammenhängen, gekoppelt sind. Ein solches System von zwei gekop- 
pelten Differentialgleichungen erster Ordnung hat die Form 


f(x,y, z) 
8(x,y,2) - 


y' 
z!’ 


Für die beiden Gleichungen werden jetzt zwei Anfangsbedingungen 
benötigt: y= Yo und z = z, für x = %. Gesucht werden die 
Kurven y(x) und z(x) „ die dem Differentialgleichungssystem 
gemigen und außerdem die beiden Anfangsbedingungen erfüllen. 
Ein System von zwei gekoppelten Differentialgleichungen erster 
Ordnung kann in einer Differentialgleichung zweiter Ordnung zu- 
sammengefaßt werden, in der neben der ersten Ableitung auch 
noch die zweite vorkommt. Der höchste Differentialquotient be- 
stimmt immer die Ordnung einer Differentialgleichung. Umgekehrt 
läßt sich eine Differentialgleichung zweiter Ordnung durch zwei 
sekoppelte Differentialgleichungen erster Ordnung ersetzen. 
Durch die einfache Substitution y' = z geht nämlich eine Dif- 
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ferentialgleichung zweiter Ordnung y" = f(x,y,y') über in die 
beiden Gleichungen 


2 
f(x,yız2) » 


y' 
zt! 


Es genügt also, das Lösungsverfahren für zwei gekoppelte Diffe- 
rentialgleichungen erster Ordnung zu kennen. 


Beschreibung des Verfahrens 


Zur Lösung der Doppelgleichung wird das Verfahren nach RUNGE- 
KUTTA jetzt aus zwei parallel verlaufenden, gekoppelten Rechmungs- 
gängen aufgebaut. Es werden also je vier Zuwächse, kı und 1, 
(I= 1 pis 4), in y- und in z-Richtung errechnet und daraus zwei 
Mittelwerte k und 1 bestimmt. Das Gleichungssystem an einer be- 
liebigen Stelle Kr In» 2 schreibt sich: 


n 
Zuwachs 
kı = 2254) 1, = %els,) 
SE 5 Ind Zu 
k, = 2-f(s,) 1, = 3 els,) 
3 ı 11.7) 
Sp: xn +75 In + k, 2, +4 . 
K, =5 f(s,) 1; -7 e(s,) 
h. . 
83: ur73 tb; Zu + 12 
h h 
ky =>f(s,) 1, =zsls,) 
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Mittelwert 


k = (K+2k,+2k,+k ) (1.8) 


3 


wi 


er I H+zly+e2l,+l,)» 


die Wertepaare für den nächsten Punkt 


Karl ” + h 
Ya MmMrK (1.9) 
Zur = Zu + 1 


dienen als Anfang für die nächste Schrittberechnung usw. 


Schrittsteuerung und Korrektur 


Auch für Systeme von Differentialgleichungen bleibt beim RUNGE- 
KUTTA-Verfahren die richtige Schrittbemessung von entscheiden- 
der Bedeutung. Der in 1.1.2 beschriebene Weg einer näherungs- 
weisen Fehlerabschätzung ist auch hier gangbar: 


Korrektur der Feinrechnung für x = x + h+h 
in y-Richtung 
(1.10) 


korrigierte Werte 


(1.11) 


7 = Ynın * 5dy 


z 


Zu4n + ö, ° 


Zur Schrittsteuerung wird der größere Betrag der beiden Korrek- 
turwerte herangezogen 


5 = max(lö,l ‚lö,|) 


und mit ihm nach (1.6) verfahren. 
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2 _%o 
Yor 20 


s„(7) 
nach (1.7) 


k(I), UT) 
nach (1.7) 


Abb. 10 Differentialgleichung zweiter Ordnung 
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Programmablaufplan 


Der Ablaufplan für zwei gekoppelte Differentialgleichungen erster 
Ordnung ist ähnlich Abb.4 aufgebaut. Die beiden Funktionen 
f(x,y,z) und g(x,y,z) werden als Unterprogramme in das Haupt- 
programm eingefüst, so daß das Hauptprogramm für Differential- 
gleichungen der unterschiedlichsten Form gilt (Abb.10). Der Ab- 
laufplan mit Korrektur und Schrittsteuerung nimmt dann zur Ver- 
einfachung hierauf Bezug (Abb.11). 


1.252 Beispiel. Einschalten eines Gleichstrom-Motors mit 
Ankerrückwirkung 


Am Beispiel 1.1.3 wurde der Einfluß der Eisensättigung auf den 
Einschaltstrom einer Spule mit Eisenkern gezeigt. Nun soll der 
Grobanlauf (Strom i und Drehzahl n) eines fremderregten Gleich- 
strom-Motors (Abb.12) gerechnet werden. Sobald durch den Anker 
der Maschine ein Strom i fließt, tritt, wegen der gekrümmten 
Eisenkennlinie, eine Schwächung des Ei Hauptfeldes 
ein, die sogenannte Ankerrückwirkung . Diese hat die Form der 
Abb.13 und kann durch die Funktion 


» 4 Schwächung 
I 


+ U = 


Abb. 12 Schaltbild eines fremderregten Gleichstrom-Motors 


1) Nürnberg, Die Ankerrückwirkung der Gleichstrom-Maschine, 
Jahrbuch der AEG Forschung 1941, Bd.8 
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-60 -40 -20 20 40 60 A 


Abb. 13 Feldschwächung durch Ankerrückwirkung 


approximiert werden. Die weitere mathematische Betrachtung führt 
auf eine "mechanische" und eine "elektrische" Differentialglei- 
chung erster Ordnung, die miteinander gekoppelt sind: 


U-Ri-cw. 


Die Berechnung von 


Stron ij A udd 
Drehzanl n= I. „u. 
7 min 


erfolgt für eine 2 kW - Maschine mit den Zahlenwerten 


Spannung U = 225 V 
Ankerwiderstand R = 3,8 Q 
Ankerinduktivität L = 0,18 H 
Trägheitsmoment I = 0,05 Ws? 


„20 U _ Ur 
Leerlaufärehzal n, = 7 ee,” 1790. 
ce = —u2 Vs 
' +(20) 
40 
Anfangswerte @® =0 und i=0 für t=0 


erste Schrittweite h 0,005 s. 
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Das Ergebnis zeigt Abb.14. Zwischen dem Anlauf mit und dem An- 
lauf ohne Ankerrückwirkung ist ein deutlicher Unterschied fest- 
stellbar. 

Gerechnet wurde mit Korrektur und Schrittsteuerung auf einem 
HP 67 nach dem folgenden Programm 


Hauptprogramm Unterprogramme 
f(x,y,z) &(x,J,2) 
mit ohne 
Ankerrückwirkung 
IBL A STO 8 ABS IBL B IBL B IBLC 
3 GSsB C x<y ns POS PO 
ST ı RCL 0 xny 1 RCLD RCL O0 
RCL 1 x EEX ; RCLC RCLC 
STOA STO 5 CHS 2 x RCL 1 
-X- GSB 0 3 RcCLC RCL 3 x 
RCL 2 GSB 0 x<y 4 : u 
STO B 2 GTO 3 [6) POS RCLD 
STO 6 STOx8 6 : RTN RCLB 
RCL 9 x 6 2 x 
x RCL 0 : * a 
—X- STO+1 x>y 1 RCL 2 
RCL 3 R GTO A + : 
stoc GSB 0 2 : PGS 
STO 7 STO-5 STO:0 sto D RTN 
-X- RCL 8 GTO A RCL C 
GTo 2 STO-4 IBL O0 x 
LBL 1 RCL 3 RCL 3 
RCL O STO:4 + : 
2 STO:5 stoc POS 
x RCL O0 RCL 8 RTN 
STO-1 STO+1 RCL 2 
STO-1 RCL 4 + 
sto 0 STO+2 STO B 
RCL 6 RCL 5 RcL 1 
RCL 2 STO+r3 RCL O 
STO 6 DSZ + 
Rt} GTO(i) sSTo A 
STO 2 RCL 2 GSB B 
RCL 7 RCL 6 RCL O 
RCL 3 STOo 2 x 
STo 7 - sto 8 
Rt 1 STO+4 
STO 3 5 STO+4 
IBL 2 : GSB C 
RCL 1 STO-2 RCLO 
STO A ABS x 
RCL 2 RCL 3 STO+5 
STOo B RCL 7 STO+5 
RCL 3 STO 3 RTN 
stoc _ IBL 3 
GSB B 4 
RCL O 5 STO:O 
x B GTOo A 
STO 4 STO-3 


1.2 Gekoppelte Differentialgleichungen erster Ordnung, Differentialgleichungen zweiter Ordnung 35 


1000 


Abb. 14 Anlaufkurven eines fremderregten Gleichstrom-Motors 


Speicherplan 
A B c D E 
Wertepaare t [9) I c frei 
1 2 3 4 5 6 7 8 9 
30 
" 


t 0 i belegt belegt belegt belegt belest 


Anfangswerte 


nl © 
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1.3 _ Gekoppelte Differentialgleichung zweiter Ordnung, 
Differentialgleichung vierter Ordnung 


1.3.1 Problemstellung und Lösungsweg 


Die mathematische Beschreibung dynamischer Vorgänge führt häufig 
auf ein System von zwei oder mehreren Differentialgleichungen 
zweiter Ordnung, die miteinander gekoppelt sind: 


f(t,y,y',2,2' ) 
e(t,y,y',2,2') » 


y" 
z" 


Für die beiden Gleichungen werden jetzt 2 x 2 = 4 Anfangsbe- 


dingungen benötigt: y=Yy y'= y, ee ZIENZL. gg, Zn Zt 


o’ o 


zur Zeit t=0,. 
Um Fragen nach einer Balkenbiegung zu beantworten, muß eine Dif- 
ferentialgleichung vierter Ordnung gelöst werden, in der Glieder 
bis zur vierten Ableitung auftreten: 

ya Eix,yıy' ,y",y" ) 
Ähnlich wie im Kapitel 1.2 führt die einfache Substitution y" = z 
eine Differentialgleichung vierter Ordnung auf zwei gekoppelte 
Differentialgleichungen zweiter Ordnung zurück: 


y" = zZ 


z" f(x,y,y',2,2') 


Beide Frobleme können also nach dem gleichen Verfahren behandelt 


werden. 


Beschreibung des Verfahrens 


Für das Rechnen mit programmierbaren Taschenrechnern ist es 
zweckmäßig, die beiden Differentialgleichungen zweiter Ordnung 
nicht weiter in vier Differentialgleichungen erster Ordnung zu 
zerlegen, sondern die Doppelgleichung nach einem von NYSTRÖM 
vorgeschlagenen RUNGE-KUTTA-NYSTRÖM-Verfahren zu behandeln. 
Dann umfaßt das Programm nur zwei parallel verlaufende, gekop- 
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pelte Rechnungsgänge mit je vier Zuwächsen, k] und 1, (I=1 
bis 4), in y- und in z-Richtung. Aus diesen Zuwächsen werden 
jetzt allerdings vier Mittelwerte k, k', 1, 1! gebildet, umy, 
y', z, und z' nach einer Schrittweite h für den nächsten Punkt 
angeben zu können. Die für Taschenrechner vorgeschlagenen Ver- 
fahren sehen sich damit von der Differentialgleichung erster 
Ordnung bis zur Differentialgleichung vierter Ordnung sehr ähn- 
lich. Sie unterscheiden sich nicht durch wachsende Schwierigkeit, 
sondern nur durch den Umfang der notwendigen numerischen Rechen- 
arbeit. 

An einer beliebigen Stelle Kr In’ Yn» Zur z, schreibt sich, 
für zwei gekoppelte Differentialgleichungen zweiter Ordnung, das 
Gleichungssystem nach RUNGE-KUTTA-NYSTRÖM 


Zuwachs 
kı = 2f(s3) 1, = 2e(s}) 
SE 
h 
k, = $2(8,) 1, =gsls,) 
k 
h 1 
Sp) + ne 0 Yn+rkıs 
1 
Zn + #z: + 5 ); z4 +1, 
h h 
= ==. 1, 
k, =$-f(s,) 1, =28(s,) (1.13) 


> 
N) 
vol 
H 
m 
u 
Dn 
m. 
+ 
P 7 
u 
IS] ey 
{ur 
ut 
a 
> 
2 


z_ + h(z} + 1,) zu + 21, ; 
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Mittelwert 
N 
u 1 2 3 
1 
Key + 2, + 2ky + ku ) (1.14) 
1 = I 1} +1y+1,) 
1 
1'=3 (15 +21, +21, +1, : 


die Wertepaare des nächsten Punktes 


mr h 

Ya mMthivy + (1.15) 
Yarı = Yn + K 

ER ea h(z! + 1) 

Zur zur 


dienen wieder als Anfang für die nächste Schrittberechnung. 


Programmablaufplan 


Das Prinzip des Ablaufplanes der Abb.10 kann für zwei gekoppelte 
Differentialgleichungen zweiter Oränung übernommen werden (Abb.15). 
Die richtige Schrittbemessung bleibt unverändert wichtig. Das 
umfangreiche Gleichungssystem wird aber für eine Fein- und eine 
Grobrechnung die beschränkte Speicherkapazität eines Taschenrech- 
ners langsam übersteigen. Auf den um Korrektur und Schrittsteu- 
erung erweiterten Ablaufplan wird daher verzichtet. An Hand der 
bisherigen Beispiele wird es sicher dem Leser im Bedarfsfall ge- 
lingen, den erweiterten Ablaufplan selber zu erstellen, indem er 
für die Schrittsteuerung den größten Betrag von jetzt vier Kor- 
rekturwerten benutzt. 


1.3.2 Beispiel. Wasserspiele 


Aus einem Wasserspeier spritzt, unter einem gegebenen Winkel « 
und mit einer bestimmten Geschwindigkeit v, ein Wasserstrahl 
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h 
or Yor Y9 
Zos Z* 


nach( 1 .13) 


Erhöhe 
I um 1 


k(I), 1(I) 
nach(1.13) 


k,k', 1,1' 


nach(1.14) 


X, Y, Y', 
z, z' 
nach( 1.15) 


Abb. 15 Differentialgleichung vierter Ordnung 
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Abb. 16 Wasserfontäne 


(Abb.16). Gesucht wird die Bahnkurve der Fontäne unter Berück- 
sichtigung des Luftwiderstandes. 

Wir wollen einen Wassertropfen als kleine Kugel ansehen und die 
Aufgabe mit folgenden Zahlen rechnen: 


Kugelradius r = 1,2.10”3 m 
2 2 2 
"Schattenfläche" A, =ır m 
Volumen V= Irr? m? 
2 
Masse m ca 
Widerstandsziffer ce = 0,09 
Ns? 
Dichte der luft P = 1,2066 IF 
Erädbeschleunigung g = 9,81 "5 
s 
Anfangswinkel @_ = 30 ° 


ulB 


Anfangsgeschwindigkeit v = 10 


Die Fragestellung führt auf zwei gekoppelte Differentialgleichun- 
gen zweiter Ordnung mit dem Ansatz der Beschleunigung in x- und 
y- Richtung (Punkte sollen Ableitungen nach der Zeit bedeuten): 
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mit dem Luftwiderstand 


4 2 . a-w 
W=zcopA,Vv = av N. 


Wir ersetzen die Bahngeschwindigkeit und den Bahnwinkel des Was- 
sertropfens durch 


v= V* + y2 und a = eroten(X) 


und erhalten abschließend 


i=- 2.0? + 32): cos [eretan(2)] 
x 
y=- 2.5? + y°)-sin[aretan(2]] -E. 


Die Anfangswerte zur Zeit t# = O sind 


= 8,6603 


{a} 
u 
< 
°, 
ir} 
°© 
ua 
R 
I 
uls ulB 


Yo v,'sin Rn 5,0000 
Schrittweite 


h 


0,01 m 


3,3936 10”? 


2a 
m 


Bl- 


Ergebnis: 


B 
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--22-200000 
-_ =. ru 
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VDODAINPW—DOD 
suse,0u0n 
PO -JDN 
WMDOPRWD NN 

oO 

“ 

o\ 

u 


42 1. Differentialgleichungen mit Anfangswerten 


Die Mittelwerte k, k' und 1, 1' wurden bei einem HP 67 parallel 
auf den Primär- und Sekundärspeichern gerechnet 


Hauptprogramm Unterpro- 
gramme für 
f und g 
LBL A SF Oo GSB B RCL 1 RCL 1 LBL B 
RCLO GSB B GSB O0 STO C + RCLD 
RCL 3 GSB O0 STO+r3 PGS RTN RCL B 
RCL 4 GSB C STO+5 GSB 3 LBL 3 : 
PS Pos GSB C POS 3 TAN! 
STOo 2 GSB O0 PG GTo a STO:4 STOo 7 
STO D GSB 1 GSB oO LBL O0 STO:5 cos 
R4 STOD STO+3 RCL O RCL 5 RCLD 
STOo 1 1 STO+5 x RCL 2 2 
STO C GSB 2 GSB 1 STO 3 + x 
Rh STO C STOD F?O STO 2 RCLB 
STO O0 POS 2 STO+4 lastx x 
ILBL a GSB 1 GSB 2 STO+5 RCL4 ei 
[6] STO B STO C RTN + RCL 6 
STO 4 1 POS LBL 1 RCL O x 
STO 5 GSB 4 GSB 1 RCL 2 x STO 
POS GSB B STO B RCL 3 2 = 
STO 4 GSB O0 2 + x RTN 
STO 5 STO+5 GSB 4 RTN STO+1 
RCLE GSB C CF o LBL 2 RTN 
-I- POS GSB B xny LBL 4 IBL C 
RCL 2 GSB O0 GSB O0 RCL 3 GSB 2 RCL 7 
STO B STO+5 GSB C 2 STOA SIN 
RCL 1 GSB 1 PGS : RCL O RCL 8 
STO A STO D GSB O0 - RCLE x 
—- POS GSB 3 x + RCL 9 
RCL C GSB 1 RCL 2 RCL O STO E - 
—- STO B STO D x RTN RTN 
Speicherplan 
A B C D E 

Wertepaare x x y y t 

0 1 2 3 4 5 6 1 8 9 

a x, x, T; Y, belegt -4 belegt belegt 9,81 

Anfangswerte 
10 11 12 13 14 15 16 bis 19 


belegt belegt belegt belest belegt belegt frei 


2. Differentialgleichungen mit Randwerten 


21 Einführung 


Bei der bisher behandelten Aufgabenstellung waren neben der Dif- 
ferentialgleichung noch Anfangswerte gegeben, d.h. gewünschte 
Daten der Lösungskurve y(x) und deren Ableitungen an einer Stel- 
le, der Anfangsstelle x, Bei vielen technisch interessanten 
Fragen müssen aber bestimmte Funktionswerte und deren Ableitun- 
gen an mehreren Stellen der Lösungskurve gefordert werden, meist 
am Anfang und am Ende eines Intervalles, d.h. an dessen Rändern. 
Jetzt sucht man also die Lösung der Differentialgleichung, die 
zusätzlich auch noch die gegebenen Bedingungen an den Rändern 
eines Intervalles erfüllt (Randwertaufgabe). So ist z.B. die 
Frage nach der Durchbiegung eines auf zwei Stützen gelagerten 
Balkens eine typische Randwertaufgabe. Die Forderung, daß an den 
Stützstellen die Durchbiegung Null sein soll, sind die Randbe- 
dingungen. 

In diesem einfachen Falle sind für die Lösungskurve y(x) die An- 
fangs- und Endordinaten y(x,) und y(x,) gegeben. Die Lösungs- 
kurve wird im Intervall 


n<r<m 


gesucht. Die vorgegebenen Randbedingungen am Intervallanfang 
werden allgemein als Anfangsrandbedingungen (ARB), die am Inter- 
vallende als Endrandbedingungen (ERB) bezeichnet (Abb.17). 
Randwertaufgaben finden sich erst ab Differentialgleichungen 
zweiter Ordnung. Außerdem sind viele von ihnen linear, d.h. so- 
wohl in der Differentialgleichung als auch in den Randbedingun- 
gen treten die Funktion y und deren Ableitungen nur in linearer 
Form auf. Linearität bedeutet aber, daß sich die Lösung aus meh- 
reren Einzellösungen zusammensetzen läßt 


y(exX) = yılz) + KoyalxX) + Kyala) + 
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ERB 


Abb. 17 Randwertaufgabe 


Dadurch eröffnet sich die Möglichkeit, lineare Randwertaufgaben 

auf Anfangswertaufgaben zurückzuführen. Da Randwertaufgaben ei- 

nen höheren Schwierigkeitsgrad als Anfangswertprobleme besitzen, 
so bescheiden wir uns für Taschenrechner auf die Lösung linearer 
Randwertaufgaben. 


2.2 Differentialgleichung zweiter Ordnung 


Eine vorliegende Differentialgleichung zweiter Ordnung 
y" = f(x,y,y') 
betrachten wir in einer etwas ausführlicheren Schreibweise: 
y"=py+tay+r ; (2.1) 
pP, 4 und r können Funktionen von x sein. Die Randbedingungen 
setzen sich im allgemeinsten Falle aus einer Linearkombination 


von Ordinaten und Kurvensteigungen zusammen: 


A 


ARB: a,y(x,) + a,y'(x,) 
(2.2) 
B . 


ERB: Bıylzy) + B,y'(x,) 
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Bei einer Differentialgleichung zweiter Ordnung beschränkt sich 
der Lösungsansatz auf zwei unabhängige Funktionen y4 und Yo: 


y(xX) = yılx) + Kyy(x) » (2.3) 


Die Bestimmung dieser beiden Funktionen über eine Aufgabenstel- 
lung mit Anfangswerten gelingt unter den folgenden Bedingungen: 


1. y, mıß die Differentialgleichung y" = f(x,y,y') erfüllen. 
2. Ya muß den homogenen Teil der Differentialgleichung 
y" = p-y' + q4-y = fuom(XrYrY') 


erfüllen (r = 0). 
3. y muß für jedes beliebige K die ARB erfüllen. 
Die dritte Bedingung ist dann eingehalten, wenn für die Teillö- 


sung yı 


(2.4) 


l 
> 


a1yı(X,) + aoyilx,) = 
und für die Teillösung Ya 
AyYo(xX,) + aoyälx,) = 0 
eilt. Der bis dahin noch freie Faktor K kann dann so gewählt 
werden, daß die Gesamtlösung y(x) auch noch die ERB erfüllt. 


Die ausführliche mathematische Beschreibung findet man in der 
Spezialliteratur, z.B. /5/. 


Beschreibung des Verfahrens 
Die Differentialgleichung zweiter Ordnung 
y" = f(x,y,y') 


wird zunächst in zwei gekoppelte Differentialgleichungen erster 
Ordnung überführt: 


y'-z 


z! pz +qy+r= f(x,y,2) » 
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Beide Gleichungen werden im Intervall %,sS X< x, zweimal mit 


der konstanten Schrittweite 


or” Ha 


er N ’ 


N ganze Zahl, 

als Anfangswertaufgabe nach RUNGE-KUTTA durchgerechnet (Abb.18). 
Erste Durchrechnung für die Teillösung y4: 

Die Differentialgleichung erhält die Anfangswerte 


y(2,) en und yilx,) = zı(%,)=0 , (2.5) 
oder, falls a, =0 ist, 
e _A 
yı(x,) = 0 und yi(x,) = 2, (x,) = % . 


Gespeichert oder notiert werden die hierbei errechneten Funktions- 
werte y,(x) sowie die erreichten Daten am Intervallende x,: 
y7(%,) und y1(x,) = z1(%,)-» 


Zweite Durchrechnung für die Teillösung Ya: 

Für diese Rechnung wird, bei gleicher Schrittweite h, nur der 
homogene Teil der Differentielgleichung f, m genommen (r=0). 
Die Anfangswerte sind 


y.(x,) =0a, und y5(x,) = 2,(x,) =-0 - (2.6) 


Y2 


Xa % 


Abb, 18 Zur Lösung einer Randwertaufgabe 
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Auch hier werden die errechneten Funktionswerte Y,(x) gespei- 
chert oder notiert, sowie die Daten am Intervallende Xy: Yo (xy) 
und y5(x,) = 25 (x,).» 

Die Anfangswerte nach (2.5) und (2.6) in der ersten und zweiten 
Durchrechnung erfüllen die Gleichungen (2.4) in einfachster Wei- 
se, wie man sich leicht überzeugen kann. 

Im anschließenden Rechnungsgang wird die Gesamtlösung an die ERB 
angepaßt, durch die Festlegung der noch freien Konstanten K. Aus 


Ke[Bayz(m) + Bo2z(m)) + Bıyılzy) + Bazıly) = 2 


_B- (Bayılp) + Bozılp) 
u FÄCHER FIAEN j 


Aus dieser Gleichung läßt sich ablesen, daß die Randwertaufgabe 
nur dann eindeutig lösbar ist, wenn der Nenner nicht zu Null wird. 
Zum Abschluß werden die notierten oder gespeicherten Funktions- 
werte zur Gesamtlösung (2.3) zusammengesetzt 


folgt 


y(x) = yılz) + Kyn(X) . 


Programmablaufplan 


Zur Vereinfachung sind im Ablaufplan Abb.19 die Berechnung der 
Mittelwerte k und 1 nach RUNGE-KUTTA als Unterprogramme gezeich- 
net. Es ist darauf zu achten, daß bei der zweiten Durchrechnung 
nur der homogene Teil der Differentialgleichung f, 


hom 
werden darf. Bei Taschenrechnern, die mit Flags ausgerüstet sind, 


genommen 
läßt sich dieses durch Setzen eines Flags bequem erreichen. Flags 


sind elektronische Schalter, die man in einem Programm entweder 
setzt (einschaltet) oder löscht (ausschaltet). 


2.3 Beispiel. Seil mit Schneelast 


Die physikalische Theorie eines durchhängenden Seiles führt auf 
die Differentialgleichung 


S-y" = -.alx) 
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Anfangs- 
werte nach 


(2.5) 


X+h, 


y+k, z+l 


speichern 


Anfangs- 
werte nach 


(2.6) 


x+h, 
y+k, z+l 
speichern 


Abb. 19 Differentialgleichung zweiter 


ja 


Ablaufplan 
--7 ähnlich 
Abb. 10 


Yılzy)» Yalx,) 
z,(x,), zo(x,) 


Ordnung mit Randwerten 
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Abb. 20 Seil mit Schneelast 


Seilzus S N, ZIast pro längeneinheit « n . 


Auf dem Seil liege eine Schneelast veränderlicher Stärke (Abb.20). 


Wir approximieren diese last durch die Funktion 


Die Berechnung des Seildurchhanges y(x) erfolgt mit dem Zahlen- 
wert 


‚ einer Schrittweite h = 0,125 m 


und den Randwerten 


ABB x,=0 m, y(x,) =0 m; also a, =1, a,=0, A=0O 
ERB x,=1 m, y(x,)=1 m; also ßı=1, ß2=0, B=1. 


Das Rechenprogramm enthält keine Ausgabe der Teillösung y, und 
Y>- In der Ergebniszusammenstellung sind sie aber mit aufgeführt. 
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x 
m 


0,000 
0,125 
0,250 
0,375 
0,500 
0,625 
0,750 
0,875 
1,000 


yı 
m 


0,000 
0,003 
0,024 
0,081 
0,190 
0,365 
0,618 
0,957 
1,388 


Programm für einen HP 


Hauptprogramm 
LBL A STOo+1 
RCL 1 Rt} 
STO A GSB 0 
1 STO-5 
- ROCL 6 
x= STO-4 
SF 2 

RCL 2 STO0:4 
STO B STO:5 
F? 1 RCL O 
GTO 3 STO+1 
GTO 4 RCL 4 
IBL 5 STO+2 
RCL 3 RCL 5 
STO C STO+3 
F? 2 GIO A 
GTo 1 LBL 0 
GSB B RCL 3 
RCLO + 

x STO C 
STO 4 RCL 6 
STO 6 RCL 2 
GSB C + 

RCL oO STO B 
x RCL 1 
STO 5 RCL O 
GSB O + 

GSB O0 STO A 
2 GSB B 
STOx6 RCL O 
x x 
RCL O STO 6 


Vor dem Programmstart ist Flag 1 zu setzen! 


Speicherpla 
Die Werte x 


n 


vo” 


= Iund x 


a0 


CF Oo 
GTO A 
IBL 2 
RCL B 
RCL 8 


x 
RCLC 
RCL 9 
x 


+ 
STO-7 
(6) 


STO 1 
RCLE 
STO 2 
RCLD 
CHS 
STO 3 
SF O 
cr 1 
GTO A 
LBL 3 
ISZ 
STO(i) 
GTO 5 
IBL 4 
F? oO 
GTO 5 
I5Z 
RCL 7 


x 
RCL(i) 


K = - 2,388 


Unterpro- 
gramme für 
f und g 


in} 
un 
Wo 

do apa ab 


sind im Hauptprogramm die Programm- 
schritte Nr.4 und Nr.88, Nr.108. 
te Steuerung ist die Zahl 9 zu speichern. 


Im Register I für die indirek- 


2.4 Differentialgleichung vierter Ordnung 5° 


[ei D E I 
Wertepaare x y z &, &n 
8) 1 2 3 4 5 6 7 8 9 
h ARB 
> #7 4 Durch: belegt belegt belegt B ß, ß, 


Bei der ersten Durchrechnung werden die Funktionswerte yı in den 
Sekundärspeichern (maximal zehn Werte) gespeichert. 


2.4 Differentialgleichung vierter Ordnung 


Das beschriebene Verfahren zur Lösung von Differentialgleichun- 
gen zweiter Ordnung mit Randwerten läßt sich auch auf Differen- 
tialgleichungen vierter Ordnung anwenden. Eine Differentialglei- 
chung vierter Ordnung 


yY = Eix,yıy',y",y") 


kann nach Kapitel 1.3 in zwei gekoppelte Differentialgleichungen 
zweiter Ordnung überführt werden. Am Anfang und am Ende des In- 
tervalls %,s x< % müssen somit bestimmte Aussagen über zwei 
ARB A, und A, und zwei ERB B} und B, vorliegen. Durch zwei 
ARB sind von den vier Funktionsdaten y, y', y", y"' am Inter- 
vallanfang nur zwei festgelegt. Demnach kann noch über die bei- 
den anderen als sogenannte freie Anfangswerte AW verfügt wer- 
den. 

In der folgerichtigen Erweiterung des beschriebenen Verfahrens 
setzt sich die Lösung der linearen Differentialgleichung vierter 
Ordnung aus drei unabhängigen Teillösungen y}, Ya und y, zu- 
sammen 


y(x) = yılX) + Koyalx) + Kyz(X) 5, (2.8) 


und zwar so, daß die Gesamtlösung y(x), für beliebige Werte von 
K, und K, ‚„ der Differentialgleichung vierter Ordnung und den 

beiden ARB genügt. Die Anpassung der Lösung y(x) an die zwei ERB 
erfolgt durch die Festlegung der bis dahin noch freien Faktoren 


K, und K, . 


52 2. Differentialgleichungen mit Randwerten 


Beschreibung des Verfahrens 
Die Differentialgleichung vierter Ordnung 


yIV = £ix,y,y',y",y") 


wird zunächst in zwei gekoppelte Differentialgleichungen zweiter 
Ordnung aufgeteilt 


y"=2z 


z" f(x,y,y',z,2') » 


Beide werden im Intervall %,s x< x dreimal mit der konstan- 
ten Schrittweite 


N ganze Zahl, 


als Anfangswertaufgabe nach RUNGE-KUTTA-NYSTRÖM durchgerechnet 
(Abb.21), 


Erste Durchrechnung für die Teillösung y}: 
Die Differentialgleichung hat die vorgegebenen Anfangswerte bei 
Nullsetzung der beiden freien AWs 


1.ARBB=A), 5 2.ARBB=A,; 1.AWV=0; 2.AN=0. (2.9) 
Notiert oder gespeichert werden die hierbei errechneten Funk- 
tionswerte y.(x) und die am Intervallende x, erreichten Daten 
b,; und bio der durch die zwei ERB gegebenen Funktionen. 


Zweite Durchrechnung für die Teillösung Ya: 

Für diese Rechnung wird, bei gleicher Schrittweite h, nur der 
homogene Teil der Differentialgleichung f om genommen. Die An- 
fangswerte sind 


1. ARB=0; 2.ARB=0; 1. AW=13; 2.AW=0,. (2.10) 
Notiert oder gespeichert werden die Funktionswerte Y5(x) und die 


am Intervallende % erreichten Daten b51 und LEP der interessie- 
renden Funktionen, 
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Y2 


o' 
IP 


Fe 01 © 1.EBB 
by, © 2.ERB 


Abb. 21 Zur Iösung einer Randwertaufgabe 


Dritte Durchrechnung für die Teillösung Y3: 
Bei gleichbleibender Schrittweite h wird nur mit dem homogenen 
Teil der Differentialgleichung gerechnet. Die Anfangswerte sind 


1.ARB=0; 2.ARB=0; 1.AN=0; 2.AW=1. (2.11) 


Die Funktionswerte y3(x%) und die interessierenden Daten by4 und 
b32 am Intervallende werden notiert oder gespeichert. 

Für die Anpassung der Lösung y(x) an die zwei ERB B, und B, 
ergeben sich zwei Gleichungen mit den beiden unbekannten Koeffi- 
zienten K und Kz: 


b241 + b31'Kz = Bi - b14 


B 


bon Ka + ben K = Bu - bin , 
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Die Auflösung dieses Gleichungssystems nach K, und K, 


(By - bandbao - (Ba - baadbaı 2 
bes (2.12) 
21 Pa2 = ba2 baı 


K, 


(Bo = Brodban = (Br = ba )bos 
Bonlb39:= Baa Pi 


i 


ermöglicht die notierten Werte der drei Funktionen Yı' Ya und Y3 
zur Gesamtlösung (2.8) zusammenzusetzen: 


y(x) = yı + Koya(xX) + Ky,(a) 


Programmablaufplan 


Im Ablaufplan (Abb.22) ist die Berechnung der Mittelwerte k, k', 
l und 1' nach RUNGE-KUTTA-NYSTRÖM als Unterprogramm gezeichnet. 
Zu beachten ist, daß die dreimalige Durchrechnung der Differen- 
tialgleichung mit jeweils anderen Anfangswerten erfolgen muß. 
Außerdem darf bei der zweiten und dritten Rechnung nur der homo- 
gene Teil der Differentialgleichung tom Verwendung finden, 

was sich durch das Setzen eines Flags wieder bequem erreichen 


läßt. 


2.5 Beispiel. Belasteter Balken auf elastischer Unterlage 


Die mathematische Beschreibung der Balkenbiegung (Abb.23) führt 
auf eine Differentialgleichung vierter Ordnung: 


Biegesteifigkeit EI, Bettungsziffer k, last pro längeneinheit q. 
Die Differentialgleichung wird überführt in 


EI: y" = z 


z"=qa-ky. 
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Abb. 23 Elastisch gebetteter Balken 


Berechnet wird die Durchbiegung eines Trägers vom Normalprofil 


NP I 20, Länge 10 m, EI = 4,5-10 


Randbe 


2. ERB 
Schrit 


dingungen 


nm, %, = 


ee: 
Mn» » 

oO 0 
De 
U [ ı 


tweite h 


Ergebnis 


oooapPnwo 3% 


Ww 
l 


y4 


‚103 
30,340 _, 
74,060 -10 


7,391: 1073 
3,993 -10*? 


1 


m, x = 1005, q = soot, 
m 
1. AW yılz,) 
2. AW 2'(x,) 
Y3 y(x) 
m mm 
0,0 +10” 0,0 
0,296 13,72 
2,370 22,00 
8,000 22,00 
18,960, 13,72 
37,030 +10 0,0 
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Programm für einen HP 67 


Hauptprogramm Unterpro- 

gramm für 
f und g 

IBL A GSB 1 GSB 4 + ISZ LBL B 

RCLO STO D cF 0 RTN cr 1 PDS 

RCL 3 1 GSB B LBL 3 GTOA RCLC 

RCL 4 GSB 2 GSB 0 3 IBL 8 RCL 6 

P9S STO C GSB C STO:4 RCLA : 

sTo 2 PS GSB 0 STO:5 STO 8 POS 

STO D GSB 1 GSB 3 RCL5 RCLC RTN 

RI STO B RCL 2 RCL 2 STO 9 

STO 1 1 STO D + 1 IBL C 

STO C GSB 4 RCL 1 STO 2 STO 4 PGZS 

Rt GSB B STO C LASTX 0) ROCL 7 

STO 0 GSB O0 PSS RCL 4 STOE RCLA 

IBL a STO+5 GSB 3 + STO 1 RCL 8 

0) GSB C FG RCL O STOo 3 x 

STO 4 GSB O0 GToO a x STO 2 F? 1 

STO 5 STO+5 IBL OÖ 2 ISZ + 

PQS GSB 1 RCL O x GTO A RTN 

STO 4 STO D x STO+1 LBL 9 

STO 5 PSS STO 3 RTN RCL A 

RCLE GSB 1 F? oO LBL 4 RCL 8 

—- STO B STO+4 GSB 2 STO:6 

1 GSB B STO+5 STO A : 

(6) GSB O0 RTN RCL O STO 5 

- STO+3 LBL 1 ROLE RCL 9 

x = STO+5 RCL 2 + x 

SF 2 GSB C RCL 3 STO E RCLC 

RCL 2 GSB O0 + RTN = 

STO B STO+3 RTN LBL 7 RCL 9 

RCL 1 STO+5 LBL 2 RCL A RCL 6 

STO A GSB 1 xny STO-6 x 

—- STO D RCL 3 RCL C STO-7 

F? 2 2 2 STO-7 R4 

GTo(i) GSB 2 : 1 STO:7 

SF Oo STO C - STO 2 RCL 7 

GSB B PQOS X je) RCL 5 

GSB 0 GSB 1 RCL O STO E x 

GSB C STO B x STO 1 STO+6 

GSB O0 2 RCL 1 STO 3 RTN 


Die zweite Durchrechnung der Differentialgleichung beginnt mit 
dem Programmschritt Nr. 147, die dritte Durchrechnung mit Nr. 161 
und die Lösung der zwei Gleichungen mit zwei Unbekannten bei 

Nr. 170. Der errechnete Wert K, befindet sich anschließend im 
Datenspeicher 6, minus Kz im Datenspeicher 7. 

Vor dem Programnmstart ist Flag 1 zu setzen! 


Speicherplan 

Die interessierenden Daten am Intervallende x“ werden im Haupt- 
programm abgefragt 

1. ERB : Nr. 143, 157, 171, 

2. ERB : Nr. 145, 159, 178. 
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Die Anfangsbedingungen für die zweite und dritte Durchrechnung 
sind im Program die Nr. 147 vis 152 und Nr. 161 bis 167. 

Der Wert = 10 ist der Schritt Nr. 22/23. Im Register I für 
die indirekte Steuerung ist die Zahl 7 zu speichern. 


A B E I 
Wertepaare y y' zZ z' x 7 
() 1 2 3 4 5 6 7 8 9 
3 ya) st(x,) za.) z’lm,) belegt 35 oO B, bg bon 
Anfangswerte E, -K, 
10 11 12 13 14 15 16 17 18 19 


belegt belegt belegt belegt belegt belegt EI q -k frei 


3. Differentialgleichungen mit Eigenwerten 


31 Einführung 


Wichtige technische Anwendungsgebiete werden durch homogene Rand- 
wertaufgaben beschrieben, Das sind Aufgaben, für die sowohl die 
Differentialgleichung als auch die Randbedingungen homogen aus- 
fallen, d.h. in denen nur Glieder mit der Funktion y oder deren 
Ableitungen vorkommen. Solche Randwertaufgaben sind häufig Eigen- 
wertprobleme, Darunter versteht man eine homogene Differential- 
gleichung, in der ein zunächst unbestimmter Parameter X auftritt, 
Z«B. 


y""-kray . 
Nur für ganz bestimmte Werte dieses Parameters X , der auch in 
den Randbedingungen enthalten sein kenn, wird die Randwertauf- 


gabe lösbar. Eigenschwingungszahlen schwingungsfähiger Systeme 
oder das Knicken eines Stabes sind solche Eigenwertprobleme. Hier 


y 


Abb. 24 Regula falsi 
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ist das Ziel der Rechnung das Auffinden der speziellen X - Werte, 
der sogenannten Eigenwerte. Damit sind die Eigenwertaufgaben eng 
mit den Randwertaufgaben des vorigen Kapitels verbunden. 

Das systematische Herantasten an den unbekannten Parameterwert, 
der eine Lösung garantiert, gelingt sehr gut mit Hilfe der ein- 
fachen Regula falsi. Diese alte numerische Methode sorst für die 
angenäherte Bestimmung der Nulldurchgänge (Wurzeln) beliebiger 
Funktionen. Zwischen zwei "falschen Ansätzen" (Schätzwerten) x, 
und x, wird die Funktion y(x) durch eine Gerade ersetzt (Abb.24). 
Der Schnittpunkt dieser Geraden mit der x-Achse wird als verbes- 
serter Schätzwert X, genommen. Auf diese Weise tastet man sich 
numerisch immer näher an den exakten Nulldurchgsang x, heran. 


3.2 Differentialgleichung zweiter Ordnung 


Die Randwertaufgabe des Kapitels 2.2 vereinfacht sich für den 
homogenen Fall: 


Differentialgleichung 

y" = f{x,y,y')=p-y' +qy (3.1) 
Randbedingungen 
ARB: a1y(x,) + ay'(x,) =0 

(3.2) 

ERB: Bıy(x,) + B,3'(x,) = 0 . 
Der Lösungsansatz reduziert sich auf 

y(x) = Ky,(x) , (3.3) 


d.h. die Eigenlösung (Eigenfunktion, Eigenvektor) y(x) läßt 
sich nur bis auf einen Faktor K bestimmen. 


Beschreibung des Verfahrens 


Die homogene Differentialgleichung zweiter Ordnung wird in der 
bekannten Weise in zwei gekoppelte Differentialgleichungen erster 
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Ordnung überführt. Dann werden beide im Intervall %,s xs xy 
mit der konstanten Schrittweite 


oh *%a N 


N ganze Zahl, 


und einem geschätzten X - Wert als Anfangsaufgabe nach RUNGE- 
KUTTA durchgerechnet (Abb.25), mit folgenden Anfangswerten: 


Ya(X,) =a, und ylx,)= zulx,)=-0, » (3.4) 


Von Interesse sind nur die mit dem A- Wert erreichten Daten 
Ya) oder ys(x%,) = 2,(%,) am Intervallende x, bzw. der Wert 
der entsprechenden Linearkombination 


Bıya(&y) + Bazp (X) =B -» (3.5) 


Bei der richtigen i/ahl von A wird B nach Voraussetzung (3.2) zu 
Null. Dieser Nulldurchgang der Linearkombination B(A) läßt sich 
mit der Regula falsi finden (Abb.26). Sie verläuft nach dem fol- 
genden Rechnungsgang: 

Aus der Differenz zweier aufeinanderfolgender Werte B, 


A=B-Bi > (3.6) 


errechnet sich der Zuwachs des Parameters X 


AN: 


(A Mar "A 9 m 


(3.7) 


Y2 


a % 


Abb. 25 Zur Lösung einer Eigenwertaufgabe 
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Br 


A 
Abb. 26 Zur Lösung einer Eigenwertaufgabe 
und damit der verbesserte X - Wert 
Ant = An + (Ada usw. (3.8) 


Der Startwert für A und AX wird geschätzt, sowie der erste, in 
der Rechnung benötigte Wert B für n = 0, B, = 0 gesetzt. Die Rech- 
nung muß durch die Vorgabe einer Genauigkeitsschranke für den Zu- 
wachs AX begrenzt werden, z.B. auf vier Stellen nach äem Komma. 
Der so ermittelte Eigenwert X ist aber nur eine mögliche Lösung. 
Weitere X- Werte findet man mit anderen Startwerten. 

Die Konvergenz der Regula falsi bei Eigenwertaufgaben ist sehr 
gut. Trotzdem kann eine Kontrollausgabe für A oder AA im Pro- 
gramm vorteilhaft sein. 


Programmablaufplan 


An einen für den homogenen Fall vereinfachten Ablaufplan der 
Randwertaufgabe Abb.19 wird die Regula falsi angehängt und die 
Rechnung mit einer Genauigkeitsschranke abgeschlossen (Abb.27). 


3.3 Beispiel. Schwingende Saite veränderlicher Stärke 


Eine elastische Saite sei an den Enden eingespannt. Sie wird an- 
gezupft. Gesucht ist die Frequenz der Schwingung. 
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Es ‚\,A A 
ı | ARB, ERB 
Anfangs- 
Start Eingaben werte nach 
(3.4) 


x+h k,l nach Ablaufplan 

k - -Jähnlich 
y+ (1.8) Abb. 10 
z+]1 


n+1 
speichern 


Adnrı nach 
(3.7) und Genauig- 
speichern r-]keits- 

I schranke 

| 

I 

nein (@) 

speichern 


Abb. 27 Differentialgleichung zweiter Ordnung mit Eigenwerten 
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Die zu lösende Differentialgleichung ist 


S.y" + X-m(x)'y= 0 


Parameter X =w® ‚ Seilzug S N, Winkelgeschwindiskeit w ze, 
Ns® 

Masse pro längeneinheit ım 7. 
m 


Die Saite soll in der Mitte doppelt so stark als an den Rändern 
sein. Wir approximieren daher die Masse m durch die Funktion 


= ; ı 
m = m,(1 + nz) ö 


Dann lautet die Differentialgleichung 


m 
ERENE, GEBEL © IE N x 
y"=- A 3 (1 + eur) s 


Die Berechnung des Parameters X erfolgt mit den Zahlenwerten 


Zo 26 + 

3 = 3,2410 5, %9,=-09 nn, %,= 1 m, Schrittweite 
m 

h= 0,1 m, den Randwerten 


ARB x 


u 
u 


a0 m ylx,)=0 m;elso au=1,0,=0, 


ERB x 


v„'T m yz)=0 m;also f,=1, B,=0, 


erster Schätzwert X = 107 ‚AA = 10° ,„ Genauigkeitsschranke 
eine Stelle nach dem Komma. 


Ergebnis: 
Parameter A= 7 242 293 = w 


Frequenz f= 37° 428 Hz 


Das Ergebnis wurde mit dem folgenden Programm auf einem HP 67 


gerechnet: 

Hauptprogramm Unterpro- 
gramme für 
f, gundB 

IBL A 2 RCL 6 RCL 9 LBL B 

RCL 1 STOx6 RCL 2 xay RCL C 


STO A x + STO 9 RTN 
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1 RCL Oo STO B - LBL C 
- STO+1 RCL 1 ST0:8 RAD 
x=0 RB RCL 0 RCL 8 RCLA 
SF 2 GSB O0 + STO+7 " 

RCL 2 STO-5 STO A RND x 

STO B RCL 6 GSB B xfo SIN 
RCL 3 STO-4 RCL 0 GTOo 2 1 

STO C 3 x RCL 7 + 

F? 2 STO:4 STO 6 RTN POS 
GTO 1 STO:5 STO+4 LBL 2 RCL 1 
GSB B RCL O STO+4 {6} PGS 
RCL O STO+1 GSB C STO 1 x 

x RCL 4 RCL O0 RCLD RCL B 
STO 4 STO+2 x STO 2 x 
STO 6 RCL 5 STO+5 RCLE RCL 7 
GSBC STO+3 STO+5 STO 3 x 
RCL 0 GTO A RTN GTOA CHS 

x LBL O0 IBL 1 RTN 
STO 5 RCL 3 GSB a IBL a 
GSB 0 + PAUSE RCL B 
GSB 0 STO C STOx8 RTN 


Die Regule falsi beginnt mit dem Programmschritt Nr. 69. Die 
Linearkombination B wird im Unterprogramm LBL a berechnet und 
zur Konvergenzkontrolle angezeigt. Die Genauigkeitsschranke ent- 
spricht dem eingestellten Anzeigeformat. 


Speicherplan 
A B C D E 
Wertepaare x y 2 an - a] 
je) 1 2 3 4 5 6 7 8 9 
$ % 9% -@, belegt belegt belegt X AA 
10 11 12 bis 19 
Mo 
frei 5 frei 


3.4  Differentialgleichung vierter Ordnung 


Für eine Differentialgleichung vierter Ordnung sind im homogenen 
Falle beide ARB und beide ERB Null: 


1. ARB=0; 2.ARB=0; 1.EB=0; 2.EB=0. 
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Der Lösungsansatz (2.8) reduziert sich auf 


y(x) = Ky-yz(xX) + Kr-yz(X) » (3.9) 


Beschreibung des Verfahrens 


Die homogene Differentialgleichung vierter Ordnung wird in der 
bekannten Weise in zwei gekoppelte Differentialgleichungen zwei- 
ter Ordnung überführt. Beide Gleichungen werden dann im Intervall 
Da < % zweimal mit der konstanten Schrittweite 


= A N ganze Zahl, 


und einem geschätzten X- Wert als Anfangswertaufgabe nach 
RUNGE-KUTTA-NYSTRÖM durchgerechnet. 


Die erste Durchrechnung für die Teillösung Ya erfolgt mit den 
Anfangswerten 


1. ARBB=0; 2. ARBB=0; 1.AW=1; 2.AN=0. (3.10) 
Von Interesse sind mur die mit dem geschätzten A- Wert am In- 


tervallende x, erreichten Daten LIE und b der Funktio- 
nen, die durch die ERB gegeben sind. 


22 


Die zweite Durchrechnung für die Teillösung Y3 erfolgt mit den 
Anfangswerten 


1.AB=0; 2.AB=0; 1.AW=0; 2. AW= 1. (3.11) 


Die am Intervallende interessierenden Daten seien b31 und b32 . 
Das Gleichungssystem zur Anpassung an die Forderung, beide ERB 
gleich Null, hat die beiden Unbekannten K, und K, 5 


b24'K, + b31°Kz {6} 


il 
(0) 


b55 Ky + b22 'Kz 
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Abb.15 


t 


speichern 


werte nach 
speichern 


Anfangs- 


Um triviale Lösungen (K, = K, = 0) auszuschließen, muß die Haupt- 


determinante zu Null werden 
D = Sbsrba - bebg=0 » (3.13) 


Dieser Nulldurchgang der Hauptdeterminante D(X) kann wieder 
mit der Regula falsi gefunden werden: 


Aus der Differenz zweier aufeinanderfolgender Determinantenwerte 
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n 


D.+1 nach 
(3.13) und 
speichern 


Rückruf D 


AAn,1 nach 
(3.15) und 
speichern 
An+1 nach 
(3.16) und 
speichern 


A 


Genauig 
keits- 
schranke 
Abb. 28 Differentialgleichung vierter Ordnung mit Eigenwerten 


errechnet sich der Zuwachs für den Parameter X 


AA 


n 
(A) 41 FR Dar (3.15) 

und damit der verbesserte XA- Wert 
Anl = An + (Ad ar usw. (3.16) 


Der Startwert für X und AX wird geschätzt, der erste benötigte 
Wert D für n=0, D = 0 gesetzt. Die Rechnung wird durch 
eine Genauigkeitsschranke für den Zuwachs AX begrenzt. Weitere 
A - Werte, falls vorhanden, findet man durch die Vorgabe anderer 
Startwerte. 
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Programmablaufplan 


An einen für den homogenen Fall vereinfachten Ablaufplan der 
früheren Randwertaufgabe Abb. 22 wird die Regula falsi ange- 
hängt und die Rechnung mit einer Genauigkeitsschranke abge- 
schlossen (Abb. 28). 


3.5 Beispiel. Knickung eines Stabes unter seinem Eigengewicht 


Wir betrachten ein Stabilitätsproblem. Ein senkrechter, unten 
eingespannter Stab konstanter Biegesteifigkeit EI werde in der 
Längsachse durch eine Kraft F belastet (Abb. 29). Die Suche nach 
der kritischen last aus Kraft und Eigengewicht, bei der eine 
Knickung des Stabes eintritt, führt auf die homogene Differen- 
tialgleichung 


EI-yN - [6,12 e x).y]' - (F + 6,1). y" 
mit dem Eigengewicht pro längeneinheit des Stabes G, Er 
Wird die äußere Kraft F = O0 gesetzt, so geht die Gleichung nach 
kurzer Umforming über in 

yY=-ıeyr ey), 
mit dem Parameter \ = 
keit. 


G 
= aus Eigengewicht und Biegesteifig- 
Dieses ist die mathematisch-physikalische Formulierung der Spruch- 
weisheit, daß die Bäume nicht in den Himmel wachsen. Die Grenze 
des Wachstums gibt die Lösung der Eigenwertaufgabe an. Sie er- 
folgt für die Zanlenwerte x =0O nn, = 1 m, Schrittweite 


a 
h= 0,1 m, den Randbedingungen 


1. ARB y(x,) =0 1. AW y'(x,) 
2. ARB z(x,) = 0 2. AW z'(x,) 
1. ERB y'(x,) = 0 
2. ERB z'(x,) = 0 


erster Schätzwert X = 10 und A/A= 1. 


Ergebnis X = 7,8377 . Weitere Durchrechnung mit anderen Stab- 
längen 1= x führt zu der Gleichung 
61? 


[e) 
Er = 1,8377 . 
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Abb. 29 
Knickung eines Stabes 


Gerechnet wurde mit dem folgenden Frogram auf einem HP 67: 


Hauptprogramm Unterpro- 

gramme für 
fundg 

IBL A GSB 2 GSB O0 IBL 3 RCL 9 IBL B 

RCLO STO C GSB C 3 x5Yy RCLC 

RCL 3 PGS GSB O0 STO:4 STO 9 RTN 

RCL 4 GSB 1 GSB 3 STO:5 - 

FQS STO B RCL 2 RCL 5 STO:8 IBL C 

STo 2 1 STOD RCL 2 RCL 8 RCLE 

STO D GSB 4 RCL 1 + POS RCLC 

R GSB B STO C STO 2 STO+7 x 

STO 1 GSB O0 Ss IASTX RND RCL B 

STO C STO+5 GSB 3 RCL 4 x£o0 + 

Ri GSBC Pg + GTOD RCL7 

sSTo oO GSB O GTo a RCL oO RCL 7 x 

IBL a STO+5 IBL O0 x RTN CHS 

(6) GSB 1 RCL © 2 IBL D PG 

STO 4 STO D x x O RTN 

STO 5 PQS STO 3 STO+1 STO E 

POS GSB 1 F? oO RTN STO 1 

STO 4 STO B STO+4 IBL 4 STO 3 

STO 5 GSB B STO+5 GSB 2 STO 4 

RCLE GSB O0 RTN STO A 1 

RCL 6 STO+3 LBL 1 RCL oO STO 2 

_ STO+5 RCL 2 RCLE 2 

x=0 GSB C RCL 3 + ST I 

SF 2 GSB 0 + STO E GTOA 

RCL 2 STO+r3 RTN RTN LBL 6 

STO B STO+5 IBL 2 IBL 5 RCLB 

ROL 1 GSB 1 xQOYy Pos STO 6 

STO A STO D RCL 3 DSZ RCLD 

F? 2 2 2 GTO 6 STO 7 

GTO 5 GSB 2 R RCLD PGS 

SF O STO C - RCL 6 je) 

GSB B POS x x STO E 

GSB O0 GSB 1 RCLO RCLB STO 1 

GSB C STO B x RCL 7 STO 3 

GSB 0 2 RCL 1 x STO 2 

GSB 1 GSB 4 + - 1 

STOD CF Oo RTN —- STO 4 

1 GSB B STOx8 GTO A 
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Die erste Durchrechnung der Differentialgleichung beginnt beim 
Programmschritt Nr. 165 (LBL D), die zweite bei Nr. 181 und die 
Regula falsi bei Nr. 161. Die Determinante D wird unter IBL 5 
berechnet und zur Konvergenzkontrolle angezeigt. Die Genauig- 
keitsschranke entspricht dem eingestellten Anzeigeformat. 


Speicherplan 


Die interessierenden Daten am Intervallende x werden im Haupt- 
programm abgefragt 


1. ERB: Nr. 146 und Nr. 177 
2. ERB: Nr. 143 und Nr. 179. 


Die Anfangsbedingungen für die erste Durchrechnung sind im Pro- 
gramm die Nr. 165 bis 174 „ für die zweite Durchrechnung die 

Nr. 182 vis 188. Das Programm sollte bei IBLD (Nr. 165) gestar- 
tet werden. Im Register I für die indirekte Steuerung ist die 
zahl 2 zu speichern. 


A B C D E I 
Wertepaare y y' 2 z'’ x 2 
6) 1 2 3 4 5 6 7 8 9 
a y(x,) y'(x,) z(x,) z'(x,) belegt x, A frei frei 
10 11 12 13 14 15 16 17 18 19 


belegt belegt belegt belegt belegt belegt b51 b AA 


22 Dn 


Teil B 


Partielle Differentialgleichungen 


4. Einführung, Gitterpunktmethode 


Fast noch zahlreicher als für gewöhnliche Differentialgleichun- 
gen sind die Anwendungsgebiete der partiellen Differentialglei- 
chungen. Durch die Möglichkeit, sie auf Rechenautomaten zu lösen, 
gewinnen sie in den letzten Jahren zunehmend an Bedeutung für 
alle Ingenieurwissenschaften. Auch programmierbare Taschenrech- 
ner sind zur Lösung geeignet. Die folgenden Kapitel sollen daher 
den partiellen Differentialgleichungen zweiter Ordnung gehören, 
Ihre Typenbezeichnungen erhalten sie nach wichtigen technisch- 
physikalischen Anwendungsgebieten: 


2 
1. . u gu 
Wärmeleitung = 2a: 
2 ze 227 
2 2. 
Wellengleichung 2 = p° . ce 
dx dt 
Telegraphengleichung eu 2, Ur 2u 
(Wanderwelle) dx° 9% d+2 
Stromverdrängung 9®u ri 9°u EEE A 
(Skineffekt) dx? dy? öt 
2 2 
: u 2 ou 
Membranschwingung <z + = b 
Ox dy dt 
Poissonsche - , 9®u 3°u fl 
: } 2° DrE xy) 
Potentialgleichung öx Iy 


Bei zwei unabhängigen Veränderlichen x und y bzw. x und t wird 
als Lösung eine Funktion u gesucht, eine Fläche im x t u- Koordi- 
natensystem, die die Differentialgleichung erfüllt und zusätzlich 
vorgegebene Anfangs- und Randwerte einhält (Abb.30). Bei drei un- 
abhängigen Veränderlichen muß die gesuchte Lösung einen Raum im 
x y t u- Koordinatensystem umfassen. Was Wunder, daß die Theorie 


76 4. Einführung, Gitterpunktmethode 


x 


Abb. 30 Fläche als Lösung einer partiellen Differentialgleichung 


der partiellen Differentialgleichungen zu den schwierigen Kapi- 
teln der Mathematik gerechnet wird. Vor diesem Hintergrunde hebt 
sich die Gitterpunktmethode (Differenzenverfahren) als eine ein- 
fache, numerische Näherungslösung wohltuend ab. Das Verfahren ist 
seit langem bekannt. Es geht auf Leonhard Euler (1707 bis 1783) 
zurück und ist auch auf programmierbaren Taschenrechnern leicht 
zu realisieren. Heute zählt es zu den wirkungsvollen Lösungsme- 
thoden, nicht nur für partielle Differentialgleichungen, sondern 
auch für nichtlineare Aufgaben der mathematischen Fhysik. Die 
Grundgedanken der Gitterpunktmethode sind: 


Das Gebiet, in dem man die Lösung u sucht, wird mit einem regel- 
mäßigen Punktegitter überzogen, welches aus gleichen, meist 
rechteckigen oder quadratischen Zellen besteht (Abb. 31). 


Die vorliegende partielle Differentialgleichung wird in den 
Gitterpunkten durch eine Differenzengleichung ersetzt. 


Abb. 31 Gitterpunktmethode 
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Aus dem Differenzenansatz entsteht ein (lineares) Gleichungs- 
system, dessen Lösung näöherungsweise den jeweiligen Wert der 
Funktion u in den einzelnen Gitterpunkten angibt. Vielfach er- 
hält man eine sehr einfach aufgebaute Rekursionsformel, in der 
ein einzelner noch unbekannter u-Wert mehreren bereits bekann- 
ten gegenübersteht. Durch wiederholtes Anwenden dieser Rekur- 
sionsformel lassen sich die u-Werte in allen Gitterpunkten fin- 
den. Eine solche mechanische Wiederholung gleichartiger Rechen- 
zyklen führt aber zu einfachen und kurzen Rechenprogrammen, wie 
es Taschenrechner erfordern. 


Für die Güte der Näherung ist es wichtig, wie der Differential- 
quotient durch einen Differenzenquotienten ersetzt wird. Man 
sollte, wenn möglich, nur zentrale Differenzenformeln benutzen. 
Das sind Formeln, in denen die Tangentenneigung in einem Gitter- 
punkt P durch die Steigung der Sehne zwischen zwei benachbarten 
Gitterpunkten ersetzt wird (Abb.32). Die richtige Wahl der 
Schrittweite h ist, wie bei den gewöhnlichen Differentialglei- 
chungen, von zentraler Bedeutung. Der Gitterpunkt P habe die 
Koordinaten 

x=i.h 


t=jk, 


die beiden benachbarten Punkte 


x+h= (ir) und x-h= (i-1).h 
t+k= (j+1)‘k und t-k= (j-1)’k 
und der Funktionswert im Punkte P sei u, ji 
’ 
u Tangente 
x 


(i-1)-h i-h (i+1)-h 


Abb. 32 Tangente und Sehne bei Zentralformeln 
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Dann wird für das Gitterpunktverfahren ersetzt 
die erste Ableitung 


du _ Mir u 
xa.n h 

du:  Misger Ya 
) = 


mit einem Fehler von der Ordnung n® bzw. x° 5 


die zweite Ableitung 


a 
dx n® 
32u a U 3410 2% j + %ı,j-1 
d+2 K 
mit einem Fehler von der Ordnung n® bzw. x? ‚ 
eine gemischte Ableitung 
au _ Mietsier T Yang T Warı,jı + Wi,jet 


9x dt 4-hk 
mit einem Fehler von der Ordnung (n+k)? : 


Schwierigkeiten können bei der mathematischen Formulierung der 
Anfangs- und Randbedingungen auftreten, weil die Zentralformeln 
an den Rändern u-Werte "jenseits" des Anfangs oder des Randes 
fordern. Sind diese u-Werte nicht bekannt, so ersetze man den 
Differentialquotienten durch eine "einseitige" Differenzenfor- 


mel 


du _ Hirt, Mi 
%- h 
du _ Mid ii 
9 k . 


Dieses ist z.B. der Fall bei der Differentialgleichung vom Typ 
"Wärmeleitung" und "Stromverdrängung". 


5. Wärmeleitung 
5.1 Problemstellung und Lösungsweg 


Die Betrachtung einer gedämpften, zeitlichen Energieströmung in 
nur einer Richtung führt auf die Differentialgleichung 


I... . (5.1) 


Ihre wichtige technische Anwendung ist die Wärmeleitung. Hier 
wird u zur Temperaturerhöhung, die vom Ort x und der Zeit t ab- 
hängt. a ist eine Materialkonstante aus Materialdichte po „ spezi- 
fische Wärmekapazität c und Wärmeleitvermögen X: 


_ eo 
wer 


Beschreibung des Verfahrens 


Die Differentialgleichung wird zunächst durch die Substitution 


t=a-r (5.2) 
vereinfacht: 

du au 

ax = ar . (5.3) 


Dann ersetzt man, wie im Kapitel 4 beschrieben, die beiden Dif- 
ferentialquotienten durch Differenzenguotienten, Meist ist aber 
die Temperaturverteilung mur zu Beginn (t = 0) der Wärmeströmung 
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bekannt, nicht aber auch noch davor. So ersetzen wir Iu durch 


oT 


eine "einseitige" Differenzenformel 


ES Ir Banken 36 Bine 55 HR ee 
2 > K 


x=ih 


r=j-k . 


Wählt man die Abmessungen der Gitterzelle h und k nicht unab- 
hängig voneinander, sondern setzt 


kenh 2 (5.4) 
so ergibt sich in den Gitterpunkten die einfache Beziehung 


- 2 u.) +Uu. . (5.5) 


a i-1,j 1,) 1,) 


i+1,j le 


%i,j+1 
Diese Rekursionsformel ist eine Gleichung zwischen der unbekann- 
ten Temperaturerhöhung u, ,j+1 im Gitterpunkt i der Zeitzeile j+1 
und drei bekannten Temperaturerhöhungen aus der Zeitzeile j. Die 
Abb. 33 soll die Rekursionsformel veranschaulichen. Der unbe- 
kannte u-Wert ist durch einen leeren Kreis, die bekannten sind 
durch volle Kreise gekennzeichnet. Die wiederholte Anwendung der 
Rekursion ergibt die Temperaturerhöhung in allen Gitterpunkten 
der Zeitzeile j+1 im betrachteten Gebiet von i=0 pbisi=1; 


Abb. 33 Rekursion: Gleichung (5.5) 


5.1 Problemstellung und Lösungsweg 


Eingaben 
Start u 
i,0 


Randwerte 


or 


Randwerte 


u N 
1,J 


Erhöhe 
j um 1 


Abb. 34 Wärmeleitung 
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anschließenä in den Gitterpunkten der Zeitzeile j+2 usw. 
Genauere Untersuchungen zeigen aber, daß dieses einfache Verfah- 
ren nur dann numerisch stabil bleibt, wenn 


n<s 0,5 


ist. Nur dann schaukeln sich unvermeidbare Rundungsfehler mit zu- 
nehmender Rechenlänge nicht auf, sondern bleiben klein oder klin- 
gen sogar ab. Andernfalls würde die Rechnung sehr schnell zu fal- 
schen Ergebnissen führen. Wählt man n -+ ‚„ so wird der Fehler 


für die Gleichung (5.5) am kleinsten, nämlich von der Ordnung ne 


Programmablaufplan 


Der Programmablaufplan (Abb.34) hat einen denkbar einfachen Auf- 
bau. Zu Beginn der Rechnung wird die bekannte Temperaturvertei- 
Jung der Zeitzeile j = 0 eingegeben. Die Temperaturen an den 
Rändern, x= 0 undx=1L (i= O0 und i = 1), müssen für alle 
kommenden j- Zeilen angebbar sein. Falls sie nicht konstant sind, 
erfolgt ihre Berechnung in Unterprogrammen. Die Temperaturerhö- 
hungen u zweier aufeinanderfolgender j- Zeilen müssen während 
der Rechnung zwischengespeichert werden. Dabei dürfen erst die 
Daten der Zeile j+1 die Datenspeicher der Zeile j-1 überschrei- 
ben. 


5.2 Beispiel. Erwärmung der Bremstrommel eines Pkw beim 
Anhalten 


Die Erwärming einer Bremstrommel (Abb.35) soll für eine Vollbren- 
sung berechnet werden. 

Geht man davon aus, daß die an den Bremsbelägen erzeugte Wärme 
radial in die zylindrische Trommel eindringt, so gilt für die 
Temperaturerhöhung u angenähert die Differentialgleichung (5.1). 
Ferner soll die Trommel die erzeugte Wärme stets aufnehmen kön- 
nen. Dieses führt zu der Randbedingung 


x=0 u = q = q,(1 - . (5.6) 


[2 


5.2 Beispiel. Erwärmung der Bremstrommel eines Pkw beim Anhalten 


Abb. 35 Bremstrommel 
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q (w/cm®) ist die Leistungsdichte an den Bremsbelägen, die, wie 
die Geschwindigkeit des Wagens, während der Bremsdauer T linear 
mit der Zeit t abnimmt. Außen, im radialen Abstand x = L, sei 


die Temperaturerhöhung auf Null abgesunken: 


gesetzt. 
Die Rechnung erfolgt mit den Zahlenwerten 


Gewichtskraft des PkW G = 14 568 N 
Geschwindigkeit v= 36 = 
Bremsverzögerung b= 9 , 
s 
Bremstrommel: 
Innendurchmesser D= 28 cm 
Breite d= 3 cm 
Materialkonstante a = 4,9 = 
cm 
Wärmeleitwert A= 0,67 ur 
cm C 


Zunächst errechnen wir die Bremsdauer bis zum Stillstand 


= on = 4,0 s 


(5.7) 


(5.8) 
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und den Bremsweg 
Ss = BLNEFT = 72 m. 
2 
Die kinetische Energie des Wagens ist 


Inn -— 1.143568 „242 _ 
W= mV =5 9,81 36° = 962 290 Jg 


und die auf ihn wirkende Bremskraft 


PH _ 2362 290 _ 13 365 N. 


Damit beträgt die Bremsleistung zu Beginn der Bremsung für alle 
vier Bremstrommeln 


P,=Fv= 13 365. 36 = 481 15 MW, 


und so erhalten wir bei einer Bremstrommelfläche von 


A= nD-d= 263,9 cm 


die anfängliche Leistungsdichte pro Bremstrommel 


_ 481 145 _ W 
Io "4.203,97 58° 


In x- Richtung wählen wir ein Gitter mit acht Punkten, bei einer 
Schrittweite h = 0,35 cn; fürn = > wird dann nach (5.4) 


k= 0,0204 cm . 


In der Randbedingung (5.6) wird die erste Ableitung durch eine 


Zentralformel ersetzt und die Substitution t=a-r (r= j-k) 
berücksichtigt 


Mo+1,j u En DR ER 2004 „ai Kk) 
= \ u . 


5.2 Beispiel. Erwärmung der Bremstrommel eines Pkw beim Anhalten = 


Nach W_4 j aufgelöst, schreibt sich schließlich die Randbedir- 
y 
sung 
2n.q . 
+ —t(1-SE;) 


%o-1,j = Mo+1,j A 7 


+ 476,2-.(1- 2) . 


= Mor1,j # 


j3 = 40 entspricht der Bremsdauer T=4 s. 
Die zweite Randbedingung für i = 1 und die Anfangsbedingungen 
für t = 0 sind nach (5.7) und (5.8) u=0. 


Abb. 36 Erwärmung einer Bremstrommel 


Das Ergebnis an der Stelle x = O zeigt die Abb. 36. Hier, an 
den Bremsbelägen, erreicht die Temperaturerhöhung den Höchst- 
wert von 330 °C nach der halben Bremszeit, d.h. nach zwei Sekun- 
den! Am Ende der Bremsperiode ist die Erhöhung wieder auf 245 x 
abgesunken, 

Der Temperaturverlauf wurde mit dem folgenden Programm auf einem 
HP 67 gerechnet, wobei sich die Ausgabe u auf die Stelle x = O 
beschränkte. 
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Hauptprogramm Ausgabe Unterprogramme 
f. die Ränder 
x=0 und x=L 


LBL A RCL(i) IBL 2 IBL a 
RCLE - RC I 1 

1 RCL(i) 1 RCLE 
+ z x£y 1 

STO E 6 RTN - 
PAUSE : RCL(i) RCL B 
16) RCL(i) -—- 5 

ST I + RTN - 
GsBa POS RCLA 
STO 0 STO(i) x 

GSB b GSsB 2 RCL 2 
STO 9 rer: + 
LBL 1 RC I RTN 
RCL(i) [e) 

1ISZ xfy IBL b 
IS2Z GTO 1 (6) 
RCL(i) PQS RUN 
+ GTO A 

DSZ 


Das Voranschreiten von j zu j+1 in den Datenspeichern erfolgt 
auf einfache Weise durch den jeweiligen Austausch von Primär- 
und Sekundärspeichern: PGS ! 


Speicherplan 
A B C D E I 
476,2 40 frei frei Zähler j zähler i 
Beginn: O 
=0 x=L 


10) 1 2 3 4 5 6 7 8 9 Zeitzeile j 
UI U Mor u u 
10 11 12 13 14 15 16 17 18 19 Zeitzeile j+1 


Die Anfangswerte für die Zeile j = O nach (5.8) sind in die Spei- 
cher 1 bis 9 einzugeben. 


6. Wellengleichung 


6.1 Problemstellung und Lösungsweg 


Ein elastischer Stab der länge L ist an seinem einen Ende fest 
eingespannt. Er wird am anderen Ende ein kleines Stück gelängt 
und losgelassen. Mit Ausnahme der Einspannstelle werden dann al- 
le Punkte Schwingungsbewegungen ausführen, die in Form einer Wel- 
le längs des Stabes ablaufen. Die Elastizitätstheorie liefert für 
die Auslenkung u eines Punktes aus seiner Ruhelage in x- Rich- 
tung (Stab-Achse) die Wellengleichung 


32u 2 au 
Ta . 
2 m. = 

Aus bb” = EA 2? 


mit Masse pro längeneinheit m ,„ Elastizitätsnmodul E und Quer- 
schnitt A errechnet sich die Wellengeschwindigkeit für elasti- 
sche Längswellen: 


Auch die Schwingung einer elastischen Saite, die zwischen zwei 
Punkten eingespannt ist, führt auf die gleiche partielle Diffe- 
rentialgleichung. Jetzt ist aber u die Auslenkung senkrecht zur 
Saite und 


mit dem Saitenzug S. 
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Beschreibung des Verfahrens 


Die beiden Differentialquotienten in (6.1) werden durch Zentral- 
formeln ersetzt: 


41 Br : SE FG el BE PER 
h k 
x=i.h 
=j k . 


Die Abmessungen der Gitterzelle h und k sollen wieder nicht un- 
abhängig voneinander sein. Die Beziehung 


n = (ui) (6.2) 


führt zu der Rekursionsformel 


= n-( + uU, (6.3) 


U; 341 wor rm je 


die aber nur für 
nsı1ı 


numerisch stabil ist. Wählt man das Gleichheitszeichen, so wird 
die Rekursionsformel besonders einfach: 


gta ig (6.4) 


Dem unbekannten Wert u, j 
’ 
Zeitzeile j und j-1 gegenüber. Die Abb. 37 soll dieses veran- 

schaulichen. 


+1 stehen drei bekannte Werte aus der 


Die Lösung u ist eine periodische Funktion. Nach einer gewissen 
Anzahl von Zeitschritten dp wiederholen sich die u-Werte, und 


die Eigenfrequenz der Schwingung läßt sich errechnen: 


1 
I = Ik Hz . (6.5) 
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Abb. 37 Rekursion: Gleichung (6.4) 


Bei einer Saite schwingen alle Punkte synchron und phasengleich, 
d.h. alle Punkte erreichen zur gleichen Zeit ihren NMaximalaus- 
schlag und gehen zur gleichen Zeit durch ihre Ruhelage. Dann dar 


gesetzt werden, und die partielle Differentialgleichung geht in 
eine gewöhnliche Differentialgleichung zweiter Ordnung mit Eigen 
werten über. 


Programmablaufplan 


Zu Beginn der Rechnung (Abb.38) werden die Anfangswerte der bei- 
den Zeitzeilen j =0 und j = 1 eingegeben. Diese können z.B. 
aus "Auslenken" und "Loslassen" bestehen. Eine vorgegebene Aus- 
lenkung liefert alle u- Werte der Zeitzeile j = 0 ; "Loslassen" 
bedeutet, daß die Rückschwinggeschwindigkeit im ersten Moment 
noch Null ist, 


oder, nach Gem Ersatz der Ableitung durch eine Zentralformel, 


,041 = 9,01. 
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Randwerte 


"0,3 


Randwerte 
ur 


Erhöhe 
ium 1 


Erhöhe 
j um 1 


Abb. 38 Wellengleichung 
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Nimmt man noch die Gleichung (6.4) hinzu, so sind die u- Werte 
der Zeitzeile j = 1 
u = I (u, + u. er 

i,1 2 ‘i+1,0 i-1,0 (6.6) 
Die Auslenkungen an den beiden Rändern x = O0 und x = L müssen 
für alle kommenden j- Zeilen angebbar sein. Die Berechnung kenn 
in Unterprogrammen erfolgen. Die Auslenkung u zweier aufeinander- 
folgender j- Zeilen müssen während der Rechnung zwischengespei- 


chert werden. Dabei darf erst die Zeile j+1 die Zeile j-1I in 
den Datenspeichern überschreiben. 


6.2 Beispiel. Träger mit Rammbär, Kompressionswelle 


Ein senkrechter Träger (Abb.39) ist am unteren Ende (x = 0) un- 
beweglich gelagert. Auf das obere Ende (x = L) fällt ein Ramm- 
bär (Masse M) mit der Geschwindigkeit 


Rammbär und Träger sollen vom Auftreffmoment an miteinander ver- 
bunden bleiben, d.h. Prellerscheinungen finden nicht statt. Für 
die Auslenkung u eines Trägerpunktes aus seiner Ruhelage gilt die 
Differentialgleichung (6.1) mit den Randbedingungen 


x=0 u=0 (6.6) 
x-L 9%u un v2. du “ b2 En M (6 7) 
z E70 o” EA 5 
und den Anfangsbedingungen 
t=0 und x Z£L u=0 (6.8) 
9u 
t=0 und x=L Te (6.9) 


Die Rechnung erfolgt für einen vier Meter langen Träger NP I 20 


Ns? 


26,81 2 


Masse pro Iängeneinheit m 


EA = 7:10 
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Abb. 39 Kompressionswelle in einem Träger 
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2 
25,48 


1,5 m . 


Rammbärmasse M 


Fallhöhe s 


Zunächst errechnen wir die Auftreffgeschwindiskeit: 


= 2:9,81-1,5 = 5,42 EE « 
Dann wählen wir ein Gitter mit acht Punkten in x- Richtung und 
einer Schrittweite h = 0,5 m. Das ergibt nach (6.2) 


26 81 es 9,8-10”? 


a ee 


In der Randbedingung (6.7) werden die Ableitungen durch Zentral- 
formeln ersetzt und k= h-b berücksichtigt: 


u 


_ Me1,j - Gl, - {?e). 2 Er nk IR a PR 
2h '$ u 


h 


Die Rekursionsformel (6.4) an der Stelle i = 1 


MEER FE Bi ES Dr BE RE a 0.10) 


liefert die zweite Gleichung für die beiden Unbekannten u, 3 
, 


und u . Nach u ausgerechnet, schreiben sich die Rand- 


1,j+1 
bedingungen 
(Arl)u 


1+1,5 
1+1,5 Z 


1-3)" 


und nach (6.6) 


nl» 
N) 
=|s 
vol 
[} 
o 
» 
a 
WW 


KR IR 6 Dale? 
"0,3 


Auch in der Anfangsbedingung (6.9) wird die Ableitung durch eine 
Zentralformel ersetzt: 


u 
1,041 161 2,2% 


2k 
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Die Rekursionsformel (6.10) an der Stelle i = 1 liefert, unter 
Beachtung, daß für alle x £{L nach (6.8) u = O ist, 


u =0 +0 - “7901 . 


Damit schreiben sich die Anfangswerte der beiden Zeitzeilen 


j=0 und ifl 4,0 * O0 mm (6.11) 
i=1 0 =0 mm 

j=1 und if£l 4,1 = O0 m (6.12) 
i=1 wı=r-r=-0,53 m. 


Die errechnete Schwingung des Trägers an der Stelle x = L und 


x= u zeigt die Abb. 39. Es ist das Bild einer Welle, die an 


den beiden Trägerenden jeweils reflektiert wird. An der Stelle 
x = L entnehmen wir als Schwingungsweite Im = 39 und finden 
damit die Frequenz nach (6.5) f = 262 Hz. 


Die Rechnung erfolgte auf einem HP 67 mit dem folgenden Programm, 
wobei die Ausgabe u sich auf die Stellen x = Y2 und x=L 
beschränkte. 


Hauptprogramm Ausgabe Unterprogramme 
f. die Ränder 
x=0 und x=L 


LBL A DSZ LBL 2 IBL a 
RCLE P93S RC I 0 

1 RCL(i) 4 RTN 

+ - : 

STO E STo(i) FRAC LBL b 
PAUSE GsB 2 xfo RCL 7 
0 PO RTN RCLA 
st ı RC I RC I x 
GSsB a 8 PAUSE RCL 8 
sTo o xfy RCL(i) + 

GSB b GTo 1 —X- RCL 8 
STo 9 PGS RTN + 
LBL 1 GTOA RCLA 
RCL(i) 2 

ISZ + 

ISZ : 


RCL(i) REN 
+ 
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Das Voranschreiten von Zeile j zu j+?! in den Datenspeichern er- 
folgt wieder durch den Austausch von Primär- und Sekundärspeicher. 


Speicherplan 
A B C D E I 
(A-1) frei frei frei zähler j Zähler i 
- 0,737 Beginn: 1 
x=0 x=L 
0 1 2 3 4 5 6 7 8 9 Zeitzeile j+1 
"0 a 


10 11 12 13 14 15 16 17 18 19 Zeitzeile j 


Die Anfangswerte für die Zeile j = O nach (6.11) sind in die 
Speicher 10 bis 18 ,„ für die Zeile j = 1 nach (6.12) in die Spei- 
cher O0 bis 8 einzugeben. 
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71 Problemstellung und Lösungsweg 


Kabel und Freileitungen der Fernmelde- und Starkstromtechnik 
sind Leitungsgebilde von großer räumlicher Ausdehnung. Es ist 
daher nicht mehr gewährleistet, daß in jedem Punkt einer länge- 
ren Leitungsverbindung zur gleichen Zeit auch der selbe Span- 
nungs- und Stromzustand herrscht, wie dieses bei konzentrierten 
Schaltungen vorausgesetzt werden kann. Jetzt het eine zeitliche 
Änderung von Spannung und Strom auch eine räumliche Änderung des 
elektrischen Zustandes auf der Leitung zur Folge. So wird die 
räumliche Verteilung des Stromes im Leitungsdraht oder die räun- 
liche Verteilung der Spannung zwischen Hin- und Rückleitung von 
Interesse, besonders bei Spannungsänderungen am Leitungsanfang. 
Die physikalische Betrachtung führt auf die beiden "Telegraphen- 
gleichungen" für die Spannung U und den Strom I: 


au _ I 
a RE = (7.1) 
en au 
= dx = G,U a (7.2) 
mit 
Widerstand der Doppelleitung R, Q/m 
Induktivität L, H/km 
Ableitung (Isolierfähigkeit) Go S/km 
Kapazität co F/ımn . 


Die beiden Gleichungen sagen aus, daß grundsätzlich jede Span- 
nungsverteilung auf der Leitung möglich ist und daß Spannungs- 
und Stromverteilung voneinander abhängen. Außerdem ist die Ver- 
teilung nicht stationär, sondern es wandern Spannungs- und Stron- 
wellen von einem Strömungspunkt aus längs der Leitung. 
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Für die weitere Betrachtung setzen wir die Ableitung G,= O0 ,„ bei 


der heutigen Isolationstechnik eine zulässige Vereinfachung. Um 
die Spannungs- und Stromwelle gleichzeitig rechnen zu können, 
führen wir eine übergeordnete Funktion u ein, aus der sich die 
Spannung U und der Strom I durch partielle Ableitungen ergeben 
sollen: 


v- u (7.3) 


x 
ou 
I=s-035 - (7.4) 


Für G,= 0 erfüllen diese Beziehungen (7.2) und überführen (7.1) 


in die partielle Differentialgleichung einer Wanderwelle: 


deu du 2 au 
— = a: + m (7.5) 
dx° x d+° 
mit 
a = RC, 
2 
bis L,C, . 


(7.5) beschreibt eine Welle, die längs einer Leitung mit der Wel- 
lengeschwindiekeit v = 1/b wandert und dabei eine Dämpfung er- 
fährt, die proportional a ist. (7.5) kombiniert die beiden par- 
tiellen Differentialgleichungen (5.1) und (6.1). 


Beschreibung des Verfahrens 


Entsprechend den Ausführungen im Kapitel 6 setzen wir für den 
Ersatz der Differentialquotienten durch Zentralformeln 


x=ih; t=jk; k=sbb=-hVic, -» (7.6) 


Dann führen wir noch zusätzlich die Abkürzung 
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ein. Dann gehen (7.3) und (7.4) über in 


a ze Se 
1, , = (w SW.) (7.8) 
is) 20 i,j-1 i,j+1 ie 
und (7.5) in 
(Artdw, 344 Et et (A-1w; ;_ı (7.9) 
mit 
L 
Wellenwiderstand 2, = Tv 2 
o 
R 
_ı_0 
und A =37 : (7.10) 
o 
In der Rekursionsformel (7.9) stehen dem unbekannten Wert w; j+ 
, 


drei bekannte Werte aus der Zeitzeile j und j-i gegenüber. 
Spannungs- und Stromwelle ergeben sich nach (7.7) und (7.8). Die 
Abb. 37 veranschaulicht auch die Rekursionsformel der gedämpften 
Welle. Die Funktionswerte in der Waagrechten ergeben die Span- 
nungswelle, die in der Senkrechten die Stromwelle. Dieser ein- 
fache Rechnungsgang erfährt leider eine Erschwernis bei der An- 
passung an die Randbedingungen, Die Verwendung von Zentralfor- 
meln führt auf zwei Gleichungen mit zwei Unbekannten. Die erste 
Gleichung ist die Rekursionsformel (7.9) für den Randpunkt, die 
zweite Gleichung ist die Randbedingung nach den Gesetzen der 
Elektrotechnik. So beschreibt z.B. die Gleichung 


R 
1 
KL a ae a Be Fe en a Pe (7.11) 
R 
1 "0,3 
ee 


Abb. 40 Speiseschaltung: Spannung mit Widerstand 


7.1 Problemstellung und Lösungsmwez =} 


nach Abb. 40 das Einschalten einer Spannung E über einen Wicder- 
stand R, am Leitungsanfang (i = 0). Die beiden Gleichungen (7.°) 


und (7.11) lassen sich entweder nach der Unbekannten wo_1 N 
’ 


(Fell I) oder nach (Fall II) auflösen. 


Wo ,j+ 
Im Fall I ist die Unbekannte ein Funktionswert jenseits des Ran- 
des (i = 0-1). Ist sie ermittelt und gespeichert, so lassen sich 
Spannung und Strom auch am Rande, wie in den anderen Gitterpunk- 
ten, nach (7.7) und (7.8) errechnen. 

Im Fall II wird die Spannung am Leitungsanfang nach (7.11) und 
der Strom nach (7.8) bestimmt. Der Rechnungsgang ist etwas um- 
ständlicher. 

Bei einem Kondensator am Anfang oder Ende der Leitung gibt es 
eine Randbedingungsgleichung in einfacher Form nur für den Fall I 
und für die Zeitzeile j+! . 


Randbedingungsgleichungen 


Eine Übersicht über einige Schaltungen an den Leitungsrändern 
zeigt Abb. 41 . Die zugehörigen Gleichungen sind für Fall I und 
II mit den folgenden Abkürzungen zusammengestellt: 


I, R} 


2 
; gs; r= (7.12) 
o ’ Bu 2, 


Schaltung a5: Einschalten einer Spannung E 


"Mo, = Mori; ” 
II (Arl)wg, jr = 2 W415 -E + (A-1)wg 5-4 
U. =E 
0,5 


Schaltung 23: Einschalten einer Spannung E über Widerstand 
und/oder Induktivität 


I (Arlırra)wo_4 5 = (Art-r-a)wo,4,; + (Art)(2aw, ; - E) + 


+ 2-(r - Ag)wo 3-1 
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R, R, 
L, L, 


Abb. 41 Randschaltungen: Speise- und Verbraucherseite 


I ENT TEN el 


E- r(w 


Or 0, a FAND je N, rg, 


Schaltung a,: Zuschalten eines geladenen Kondensators 


I (Arl)wo_4 341 Wo; + Wua,j tr A-wo_1,35-1 + 


5 Ten "or, gt ga, “| 


Die Strom- und Spannungswelle erfordert die Eingabe von w- Werten 
für die Zeitzeile j=-|1 und j=0. 


7.1 Problemstellung und Lösungsmez o 


Für 1 >0 sind die Anfangswerte w dieser Zeitzeilen Null zu set- 
zen. 
Der erste nach den obigen Gleichungen errechnete Randwert ist 


im Fall I Wo-1,0° 


im Fall II wo, 2 


Nur für Schaltung a, ist der erste Randwert w._, „ . Hier ist die 
’ 


Spannung am Kondensator zusätzlich negativ einzugeben: 


- U 


"0-1,0-1 = Wo-1,0 = "0, ° 


Schaltung b: Leitungsende offen 
w . = 2Bew. . -w r 
1+1,j 1,j-1 1-1,j 


"1,3417 #1,3-1 


U]; = 2m) 1 7 mo1,3) 


Schaltung by: Leitungsende kurzgeschlossen 


a 1,5 

II (AHldw jr = 2m 41,5 + (A-1)w] ;- 
U .= 0. 
1,5 


Schaltung b3: Leitung mit Widerstand und/oder Induktivität 
belastet 


I (Arten )wi,1 5 = 2.[atart)w ; + (Aalen, 1] i 


+ (Arl-r-a)wı_ 1; 
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II (Arlırıq)w = 2:(w + aw] 5) + (A-1+r-a)w, ;_4 


1,j+1 1-1,5 


U ee erg 


Schaltung b Leitung mit Kondensator belastet 


4° 


+ (A-1)w 


De ge FT gr 


2 _ 5 } 
en (I rn) 


Programmablaufplan 


Den Programmablaufplan für den Fall I zeigt die Abb. 42. Er ist 
wie der in Abb. 38 aufgebaut. Lediglich die Berechnung von Span- 


nung U, und Strom I; j ist hinzugekommen. 
’ , 


J 
Die Eingabe der Anfangs- und Randwerte wurde bereits bei den 
Randwertgleichungen besprochen. Die w- \lerte zweier aufeinander- 
folgender j- Zeilen müssen während der Rechnung zwischengespei- 
chert werden, Dabei dürfen die Datenspeicher der Zeile j-1 erst 
durch die Daten der Zeile j+1 überschrieben werden. 

Die Änderungen im Programmablaufplan für den Fall II sind in der 


Abb. 54 enthalten. 


7.2 Beispiel. Speisung einer am Ende offenen Leitung über eine 
Induktivität 


Es soll die Wanderwelle ermittelt werden, die beim Zuschalten ei- 
ner Spannung über eine Induktivität entsteht, wenn das Leitungs- 
ende offen ist (Abb.43). 

Die Leitung habe die folgenden Kenngrößen: 


Länge S o:= 1 km 
Be N WW 
In = 40-100 u 
R,= 16,6 Q/ım 
L_= 200-10°° H/km 


Abb. 42 


Eingaben 


Wil 


”;,c 


Rückruf 


Wanderwelle (Fall I) 


Randwerte 


wo_1 j 


Randwerte 


W141,j 


“i-1,3 


= Wir1,j 


"i,j-1 


— „speichern 


nach (7.7) 
=] und (7.8) 


Erhöhe 
j um 1 


} 
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Abb. 43 Speisung einer am Ende offenen Leitung über eine 
Induktivität 


Der relativ hohe Widerstand entspricht einem Aluminiumdraht von 
2,2 mm Leitungsdurchmesser. 

Wir unterteilen die Leitungslänge in vier Teile. Dieser sehr gro- 
be Gitterraster wurde für einen Vergleich mit Ergebnissen aus dem 
Beispiel 9.2 gewählt. Das ergibt nach (7.6) 


n=0,5 m; k=110% s. 


Die Geschwindigkeit der Welle errechnet sich zu v = 2,5-10° kn/s. 
Bei einer Länge von einem Kilometer hat die Welle das Leitungs- 
ende nach einer Zeit von 


1 _,3.40° 
T= ri 4-10 s 


erreicht. Für das Beispiel gelten die Randbedingungsgleichungen 


a7 mit r= 0 


+ (Arl)(2g-wo ; - E) + 


7 (Arla)wo_ı, 5 = (Art-a)dwo,g 5 
- EAU NG 4 
A = 0,0415 nach (7.10) 
q= 1,6 nach (7.12) 
und b} 
I. m, 29 Ben 


Die gerechnete Spannungswelle an drei Leitungspunkten zeigt die 
Abb. A4. Die Dämpfung der Welle durch den sehr großen Leitungs- 
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ohne Stromverärängung mit Stromverdrängung 


uk 


1,5 
Leitungsanfang 


1,0 


0,5 


le} 


Leitungsmitte 


1,5 


lc 


Leitungsende 
1,5 


Abb. 44 Wanderwelle des Beispiels 7.2 ohne und mit 
Stromverdrängung 
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widerstand macht sich deutlich bemerkbar. Die Spannung am Lei- 
tungsende hat einen ähnlichen Verlauf wie die Kompressionswelle 
in Abb. 39. In die Abb. 44 sind auch die Ergebnisse des Beispie- 
les 9.2 zum Vergleich eingezeichnet. 


Die Wanderwelle wurde auf einem HP 67 mit folgendem Programm ge- 
rechnet. Die Datenausgabe beschränkt sich auf Anfang, Mitte und 
Ende der Leitung. 


Hauptprogramm Ausgabe Unterprogramme 
Strom und f. die Ränder 
Spannung i=0-1; i=- 1 
LBL A RCLA IBL 2 IBL a x 
RCLE 2 RC ıI 2 POS 
1 + 2 RCLC RCL 1 
+ B : x PGS 
STO E RCL(i) FRAC RCL 1 x 
PAUSE xD x=0 x - 
(6) sTto(i) RTN RCLD RCLA 
ST I - Rl - 2 
GSB a POS x5I RCLA + 
STO 0 GSB 2 PAUSE 2 RCL C 
GSB b RC I xyI + + 
STO 6 5 RCL B x : 
LBL 1 xfy B LJASTX RTN 
RCL(i) GTO 1 —- RCL C 
152 POS ISZ - LBL b 
ISZ GTo A RCL(i) RCL 2 POS 
RCL(i) DSZ x RCL 5 
+ DSZ + RCL 5 
DSZ SPACE RCL A + 
POS RCL(i) 1 POS 
RCL(i) - + RCL 4 
RCLA IS2Z RCLC - 
x —X- x RTN 
+ RTN 2 
Speicherplan 
A B c D E I 
(A-1) 20 q E Zähler j Zähler i 
70,982 50 146 10  Begim: -1 
x=0 x=L 
0 1 2 3 4 5.6 7 8 9 Zeitzeile j 
Wo Wo Wori WW Wa frei frei frei 


10 11 12 13 14 15 16 17 18 19 Zeitzeile j+1 


8. Stromverdrängung, Skineffekt 


8.1 Problemstellung und Lösungsweg 


Eine Wärmeleitung in x- und y- Richtung führt auf die partielle 
Differentialgleichung 


au au. 
ee ee 2) 


Als Lösung wird ein Raum im x y t u - Koordinatensystem gesucht. 
Schwierigkeiten entstehen dabei weniger in der Numerik, mehr bei 
der bildlichen Darstellung der Ergebnisse. 

Neben der Wärmeleitung in dünnen Platten erlaubt die Gleichung 
die Stromverdrängung (Skineffekt) der Elektrotechnik zu berech- 
nen. Nur für Gleichspannung ist die Stromdichte über dem Leiter- 
querschnitt konstant. Bei Wechselspannung führt die endliche Auf- 
bauzeit des elektromagnetischen Feldes im Drahtinnern zur Störung 
dieser gleichmäßigen Stromverteilung. Im massiven Leiter bilden 
sich Wechselströme aus, die sich dem eigentlichen Leiterstrom 
überlagern und die gleichmäßige Verteilung aufheben. Der Strom 
wird immer mehr aus dem Innern an die Außenhaut gedrängt. Bei 
vorgegebener Spannung verkleinert sich der fließende Strom durch 
diese auch als Selbstinduktion bezeichnete Wirkung, welche einer 
Widerstandserhöhung gleichkommt. Die Funktion u ist jetzt die 
Komponente der elektrischen Feldstärke E im Drahtinnern, die pa- 
rallei zur Leiterachse gerichtet ist. Die Konstante a ist das 
Produkt aus elektrischer Leitfähigkeit K und magnetischer Per- 
meabilität u: 


a = KeN p (8.2) 
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Beschreibung des Verfahrens 


Das Verfahren aus Kapitel 5 wird in der y- Richtung erweitert. 
Zunächst vereinfacht sich die Differentialgleichung durch die 
Substitution nach (5.2) 


zu 


a2u 92u _.. gu 
FB: 7. (8.3) 


Dann werden die drei Differentialquotienten durch Differenzen- 
quotienten ersetzt, wobei man für en wieder auf eine einseiti- 
ge Differenzengleichung zurückgreift. Mit 


x=i.h 
y=zi'h (8.4) 
T= j-k= jenn® 


ergibt sich die Rekursionsformel in den Gitterpunkten 


u + (8.5) 


date" UL WEITER En 9 Po DE ee Pe 


+Uu, „, .- Au, . .) + U, , Pi 
i,i'-1,5 1,ir,3) i,i',)j 


Es stehen der unbekannten Feldstärke u, in der Zeit- 


i,i',j+1 
schicht j+1 fünf bekannte der Zeitschicht j gegenüber. Die 
Rechnung ist für ns 0,5 numerisch stabil und hat für n -+ 
den kleinsten Diskretisierungsfehler. 


Abb. 45 Zur Randwertbestimmung 


8.2 Übergang auf Polarkoordinasen 2 


Für die Bestimmung der Randwerte, der Stromdichte S unä des _ei- 
terstromes I bediene man sich der folgenden einfachen Beziehur- 
gen: 

Die Spannungsdifferenz U. zwischen zwei um d entfernte Leiter- 


punkte (Abb.45) bestimmt die Felästärke am Leiterrand 
U. 
ie 2 a u . 
Ya,j Ba; d 18:81 


Die Stromdichte S im Leiter findet man aus der Feldstärke 


S=-KU=K'E Pi (8.7) 


und den im Leiter fließenden Strom I erhält man aus der Strom- 
dichte durch eine numerische Integration über den Drahtquer- 
schnitt A (siehe hierzu Abschnitt 8.2) 


1a (ls... (8.8) 
A 


Programmablaufplan 


Der Leiter habe einen rechteckigen Querschnitt a-b „ Zu Beginn 
der Rechnung (Abb.46) erfolgt die Eingabe der Feldstärkenwerte 
der Zeitschicht j = 0 (meist gleich Null). Die Berechnung der 
Randfeldstärke Y,j aus der aufgedrückten Leiterspannung nach 


(8.6) kann in einem Unterprogramm erfolgen. Ebenfalls die Be- 
rechnung des Stromes nach (8.8). Die Feldstärkenwerte zweier auf- 
einanderfolgender Zeitschichten j müssen während der Rechnung 
zwischengespeichert werden. Die Daten der Schicht j-1 dürfen 
dabei erst durch die Daten der Schicht j+1 überschrieben werden. 


8.2 Übergang auf Polarkoordinaten 


Oft sind technische Probleme für Körper von zylindrischer Form 
zu lösen, z.B. die Stromverdrängung in runden Drähten. Zur Ver- 
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Kar) > 

s S 

- - 

= 

N = 

a - ri 
n 2 

A ; f 
7 J Pe] 
gr = u. 

5 

& 2 > 
] = B - “- 

- ee 

& La U 
Fe er m 
I - 

Ka} Ba} 


einfachung der Rechnung führt man dann gerne die Polarkoordinaten 
rund ein: 


x = r-cosy (8.9) 


y= r-sinyp . 


Nach dem Übergang auf Polarkoordinaten schreibt sich die Diffe- 


rentialgleichung (8.1) 


2 2 
deu 1, du 12.50 au 
+ "+ —- — a-77 . (8.10) 
y 
r 
) 
x 


Abb. 47 Polarkoordinaten 


8.2 Übergang auf Polarkoordinater. 


-meioer 


Abb. 46 Stromvrerdrängung 


Aus Symmetriegründen innerhalb eines Zylinders ist in vielen 
technischen Anwendungen die Funktion u unabhängig vom Winkel yp, 


d.h. 


und die Gleichung (8.10) vereinfacht sich mit der Substitution 
(5.2) zu 


au 1.9u du 
ee negr en 


Der Ersatz durch Differenzenguotienten führt mit 


(8.12) 


4 
r 
2 
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zur Rekursionsformel in den Gitterpunkten 


(8.13) 


1 
% jr ® nu, +, 43*+2 W415 - 41,5) + 
+ we 2)u 
n isj 5 
wiederum numerisch stabil für n<0,5 undn -7 als bester 


Wert, 


Jetzt werden für den unbekannten Wert u, 


i,j+1 nur noch drei be- 
’ 


kannte benötigt. Dafür muß aber gewährleistet sein, daß im Mit- 


telpunkt der Kreisfläche (i = 0) LEBE FE endlich bleibt. Die Funk- 
’ 


tion u ist symmetrisch zum Kreismittelpunkt, d.h. 

Most At, 
oder 

u - u ) 

2i 0-1,j O+1,j 
Im Kreismittelpunkt schreibt sich damit die Rekursionsformel 

1 

ug jr ® n.[2 u; en - 2uo,;] s (8.14) 
Auch die Berechnung der Stromstärke I aus der Felästärke u=E 
nach (8.8) vereinfacht sich in Polarkoordinaten. Mit dem Flächen- 


element 


dA = 2-.T-r-dr 


dr 


Abb. 48 Flächenelement in Polarkoordinaten 


8.2 Übergang auf Polarkoordinaten ı12 


y4 Ya 


Abb. 49 Keplersche Faßregel 


geht das Doppelintegral über in 


wobei a der Außenradius der Kreisfläche ist. 


Dieses Integral kann numerisch durch ein- oder mehrfachen Ansatz 
der Keplerschen Faßregel gelöst werden. Diese einfache alte Re- 


gel erlaubt ein Integral der Kurve y = f(x) angenähert aus nur 
drei Ordinatenwerten zu bestimmen (Abb.49): 


+h 
\ y-dix = 2.(y + 4: yı + Y,) 
3 (6) 1 2 5 
-h 


Wir verdanken diese Regel Johannes Kepler (1571 bis 1630). Seine 
erfolgreiche numerische Auseinandersetzung mit unlauteren Wein- 
händlern gab dieser Regel ihren Namen . . 


2) W. Gerlach und M. List, Johannes Kepler, München 1971 
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Abb, 50 Stromverdrängung, Übergang auf Polarkoordinaten 


_ 
_. 
in 


8.3 Beispiel. Stromverdrängung in runden Drähten bei nichtsinusförmiger Spanrung 


Mit y = E-r ergeben sich nach Kepler für I ,„ bei verschieder 
feiner Gitterunterteilung, folgende Gleichungen: 


Ken I-&:6,.142 
EN 

3 3 1,5 2,5 
ei Et 

ai 4-3 ] 4 
I, = Kkhy|7 (E1,; + BE, ; + Ey ,;) + Ey; (8.15) 
nee 15:0, 1,2, 94,5, 6 

5 412 Ge "E 5 ] 
I, = KA 313 (Ei; + Bj + 3 E,; + 2 E, j* 5 E, ;) + Be; 
Kreisfläche A= n.a® . 


Der Programmablaufplan der Abb. 46 vereinfacht sich zur Abb. 50. 


8.3 Beispiel. Stromverdrängung in runden Drähten bei 
nichtsinusförmiger Spannung 


Es sollen für vier Spannungskurven (Abb.51) von gleicher Frequenz 
und gleichem Waximalwert die Stromkurven und die Wechselstrom- 
widerstände eines Kupferdrahtes errechnet werden: 


Frequenz f = 13 153 Hz 
Drahtradius a = 0,33 cm 
Be aaa : 4 1 
Leitfähigkeit K= 50-10 Och 
Permeabilität = 47-107 23 


116 8. Stromverdrängung, Skineffekt 


Kurve 1 Kurve 2 


Kurve 3 Kurve 4 
Abb. 51 Spannungskurven zum Beispiel 8.3 


Der Spitzenwert der Außenfeldstärke soll jeweils 


ER 
(Edmax = 2,910 N 


betragen. Wir unterteilen den Radius in sechs Teile 
h = 0,055 em. 


Dann errechnet sich zunächst, mit r nach (8.12), der Substitution 
nach (5.2) und bei n = 1/6 


wt = 2nf-Kkun-h-j = 0,2618-j , 


d.h. die Spannungskurve hat eine Periodendauer in Zeitschritten 


: en _ 
In = a0 4 - 


8.3 Beispiel. Stromverdrängung in runden Drähten bei nichtsinusformiger Sparrurz 


Abb. 52 Stromkurven des Beispiels 8.3 
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Die Feldstärke im Drahtinnern errechnet sich nach (8.13) und 
(8.14), die sich daraus ergebenden Stromwerte nach (8.15). Das 
Ergebnis im eingeschwungenen Zustande (nach j = 120) zeigt die 
Abb. 52. Der Wechselstromwiderstand Z gegenüber dem Gleichstrom- 
widerstand R ergibt sich aus den Effektivwerten (quadratische 
Mittelwerte) von Feldstärke und Stron: 


(8.16) 


Die Stromkurven wurden auf einem HF 67 gerechnet. Im nachfolgen- 
den Programm darf jeweils nur eines der Kurvenunterprogramme be- 
nutzt werden. Die Datenausgabe beschränkt sich auf die Stromstär- 
ke I. 


Hauptprogramn Unterprogramm Unterprogramme f. Rand r-a 
I m. Ausgabe Kurve 1 Kurve 3 _ Kurve 4 

IBL A 2 IBL b IBL a IBLa IBL a 

RCLE : RCL 1 RAD RCLE ö 

1 DS2Z RCL 2 RCLE RCL B 5 

+ RC I + RCLC B RCLE 

STO E : RCL 3 x FRAC RCLB 

PAUSE - 3 SIN 4 R 

0 RCL(i) x RTN 1/x FRAC 

ST I 4 + x>y x=y 
GSB a x RCL 4 Kurve 2 GTO 3 GTo 3 

RCLA + + 3 x>y 

x 6 RCL 4 IBLa x GTo 4 

STO 6 5 + RCLE x<y 2 

GSB b POS RCL 5 RCL B GTo 4 x 

RCL oO STO(i) 5 ® xy RTN 

RCL O RC I x FRAC 4 LBL 3 

+ 5 + 1 x (6) 

RCL 1 xfy 2 0) 2 RTN 

+ GTO 1 x x - LBL 4 
3 GTO A 3 x=0 CHS 2 

s : 5 RTN x 

POS RCL 6 ENTER IBL 4 2 

STO Oo + 5 xOYy _ 

IBL 1 9 x=y 1 RTN 
Ss : 0 - 

Be RCLD x<y 4 

RCL(i x CHS x 

152 —X- 1 RTN 

ISZ RTN x LBL 3 

RCL(i) 5 xOy 

+ : 


8.3 Beispiel. Stromverdrängung in runden Drähten bei nichtsinusförmiger Spanrurg 


Das Voranschreiten in den Datenspeichern von j zu j+1 erfolgt 
durch einen Austausch der Primär- und Sekundärspeicher: FSS . 


Speicherplan 
A B c D E I 
(Bar In KA Zähler j Zähler i 


2,9-10° 24 0,2618 1,711-10° Beginn: -1 


r=0 r=a 
[6) 1 2 3 4 5 6 7 8 9 Zeitzeile j 
Eo,j Bj frei 


10 11 12 13 14 15 16 17 18 19 Zeitzeile j+1 


9. Wanderwelle mit Stromverdrängung 


91 Problemstellung und Lösungsweg 


Nach Kapitel 7 ist für die Dämpfung einer Wanderwelle der Wider- 
stand der Doppelleitung maßgebend. Im Kapitel 8 wurde gezeigt, 
daß bei schnellen zeitlichen Spannungsänderungen der Skineffekt 
eine Erhöhung des Drahtwiderstandes bewirkt. Bei der großen Ge- 
schwindigkeit der Wanderwelle ist ein Stromverdrängungseffekt auf 
der Leitung nicht auszuschließen. Seine Auswirkung auf die Welle 
soll untersucht werden. 

Setzen wir einen runden Drahtquerschnitt voraus, so ergibt die 
Zusammenstellung der Gleichungen aus Kapitel 7.1 und 8.2 ein 
System von drei partiellen Differentialgleichungen und einer In- 
tegralgleichung: 


a1 


- 5 = 25, + 1. 4 (9.1) 
91 au 

nn: = Cr (9.2) 

0®E , 1.28 _ 23E (9.3) 

dr: + r 2 T 9.3 
a 

I= 2rr-| E-r-är. (9.4) 
0 


In Gleichung (9.1) wurde der ohmsche Spannungsabfall bei Gleich- 
strom IR, durch 2-E, ersetzt, wobei En die parallel zur Lei- 


terachse gerichtete Komponente der elektrischen Feldstärke an der 
Drahtoberfläche (r = a) ist. 


9.1 Problemstellung und Losungswez er 


Beschreibung des Verfahrens 


Zur Lösung des Gleichungssystems mit der Gitterpunktmethode ver- 
wenden wir zwei Gitter: 


1. In Leitungsrichtung (Wanderwelle) 


{a} 
il 
7 
[TR 
g 


2. In radialer Richtung, senkrecht zur Leitung (Stromverdrängung) 


r= h'i' cm 


T=k'j' cm” ; k' = 4m)? R 


Dann ist das Verhältnis der beiden Zeittakte zueinander 


Ler \ (9.5) 
j' 7500-5’R, v L,C, 
Leitungslänge s km 


Widerstand der Doppelleitung R, N/ıkm 


Induktivität L, H/km 
Kapazität C, F/km 
Drahtradius a cm . 


Der Beiwert H richtet sich nach den gewählten Gitterabmessungen. 


Beiwert H 


a 
2 
2a 
ji 
a 
Ss 
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Die Berechnung der Wanderwelle soll nur während einer Zeit T nach 
dem Zuschalten einer Spannung erfolgen. Gilt für diese Zeit 


T<lLxk , 
je 


so ist der Vorgang so kurz, daß das elektromagnetische Feld nur 
wenig in den Draht eindringen kann. Die Rekursionsformel für das 
Gitter in Leitungsrichtung entspricht dann der Gleichung (7.9) 


(A'+1)w. + W224 


i,jr1 @ Tirt,j i- » (A'-1)m, 


; . . .6 
„J i,j-1 (9.6) 
A wird durch A' ersetzt, weil nicht der volle Drahtquerschnitt 
für die Stromwelle zur Verfügung steht. Je nach der gewählten 


Gitterabmessung in radialer Richtung ist 


ira At = JA 
h'=7 A' = 6A 
wer At = 9A 
zu setzen mit 
hR 
AGs R=s—g Ym. (9.7) 
° knra 


Die Wanderwelle mit Stromverdrängung wird dann weiter nach Kapi- 
tel 7 gerechnet. 
Ist dagegen 


n>_.x , 
J' 


so dringt das elektromagnetische Feld mehr und mehr in das Draht- 
innere ein. Die Welle läßt sich dann mit einem Taschenrechner auf 
folgende Weise rechnen. 


Wegen der begrenzten Anzahl von Datenspeichern verwenden wir zwei 
grobe Gitterraster 


h= km , ie cm. 


B=3 
4 


9.1 Problemstellung und Lösungswer re 


Damit sind 19 Datenspeicher belegt: je fünf für die Funktiors- 


werte w. „ und w. im Wanderwellengitter und drei al rei 


I,7 i,j+? 


für die Felästärken im Stromverdrängungsgitter Eo ‚,Eı» EM Ge 


Y 


Zeitzeile j' an den drei Leitungspunkten i = 1, 2, 3 (Abb.53). 
An den Rändern i = 0 und i = 4 werden die Feldstärken im Dreht- 
innern Null gesetzt. Diese Vereinfachung erlaubt die Benutzung 
der Randbedingungsgleichungen aus Kapitel 7. Weitere Datenspei- 
cher lassen sich bei Verwendung der Gleichungen II einsparen. 
Das Zeittaktverhältnis (9.5) wird auf die nächste ganze Zahl ge- 
rundet. 

Bei der vorgesehenen geringen radialen Unterteilung im Stron- 
verärängungsgitter sind im Drahtinnern nur zwei Feldstärken zu 
berechnen. Hierzu ersetzen wir die Gleichung (9.3) durch die Re- 
kursionsformel (8.13) bzw. (8.14). Mit u gleich der Feldstärke E 
und n = 1/6 lassen sich dann beide Felästärken für die Zeit- 
zeile j' an einem beliebigen Leitungspunkt i angeben: 


1,rr, ; 
h-E =7 (n E1 3,j'-1 + 2h Eo,i,j'-1) (9.8) 


0,i,5' 
h-E = 4. |3n-2 + 8h-E )+h-E 
1,i,j' 7 4 [3 0,i,j'-1 1,1,j'-1 a,i,j'-1) ° 


Als nächstes stellen wir die Gleichung (9.1) nach 2'E, um und 


führen die übergeordnete Funktion u „ Gleichung (7.3) und (7.4), 
ein: 


2 2 
u Tu 
2E = = + L,0 4 * (9.9) 
a ax oo dt? 
al BE ee EZ. 0,4,3° 
Stromverdrängungs- 
u Fe 
gitter u 1is) 
rza oe. oo En i,j' 
i= 0 1 2 3 4 
. eo 00 0 ee "ij 
Wanderwellengitter 


x= 0 Ss 


Abb. 53 Datenspeicher-Gitter für Wanderwelle 
mit Stromverdränguns 
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speichern 


Rückruf 


Anfanrswerte 


105 
nach (7.8) 


Nach dem Ersatz der Differentialquotienten durch Zentralformeln 
und der zusätzlichen Abkürzung 


geht (9.9) über in 


h-E tm 1 m + W Ns (9.10) 


ea a WW. . . a 
a,is) i,j+ i,j- i+r1,j i-1,j 


Anschließend wird in der Gleichung (9.4) der Strom I durch (7.8) 


ersetzt, das Integral über die Keplersche Faßregel nach (8.15) 
ausgedrückt: 


I. -2 1.02. 
ade ee 


9.1 Problemstellung und Lösurgrmer = 


[ersinern 
a Tersnern 


Abb. 54 Wanderwelle mit Stromverdrängung 


Aus den beiden Gleichungen (9.10) und (9.11) lassen sich nun an- 
geben: 

der Randwert im Stromverärängungsgitter an einem beliebigen Lei- 
tungspunkt i 

Harman) 02) 


(3Ar1)h-E, 3,5: = 31-(2-w, ı- m 


‚j- 1 


Fer dsd 
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und die Rekursionsformel für das Wanderwellengitter 


(3A+1)w, 


1,1 ® Wiri,j + "1, + (3A-1)w. - &h-E] ;,j' (9.13) 


i,j-1 


mit A nach (9.7) - 

Nun steht der Berechnung der Spannungs- und Stromwelle nach den 
Gleichungen (7.7) und (7.8) aus den Werten w ,„ wie im Kapitel 7 
beschrieben, nichts mehr im Wege. 


Programmablaufplan 


Der Programmablaufplan (Abb.54) ist eine Erweiterung der Abb.42. 
Zu Beginn der Rechnung muß das Zeittaktverhältnis (9.5) als eine 
ganze Zahl eingegeben werden. Für die Anfangswerte sei auf Kapi- 
tel 7 verwiesen. Im Rhythms des Zeittaktverhältnisses werden die 
inneren Feldstärken Eo und E, in den Leitungspunkten i=1,2, 
3 nach (9.8) neu berechnet. Der weitere Ablauf entspricht der 
Abb. 42 , mit Änderungen für die Randbedingungsgleichungen des 
Falles II. Die w- Werte zweier aufeinanderfolgender j- Zeilen 
müssen zwischengespeichert werden. Die Datenspeicher der Zeile 
j-1 dürfen erst durch die Daten der Zeile j+1 überschrieben 
werden. 


9.2 Beispiel. Einfluß der Stromverdrängung auf das Beispiel 7.2 


Das Beispiel 7.2 soll unter Berücksichtigung der Stromverdrän- 
gung gerechnet werden. Für 


h= h!'= 


rolm 


SS . 
4 


ist das Zeittaktverhältnis 


Wir benutzen ferner die Randbedingungsgleichungen a 
und b, in der Form II. 


3 nitr=0 


9.2 Beispiel. Einfluß der Stromverdrängung auf das Beispiel 7.2 127 


Die gerechnete Spannungswelle am Anfang, Mitte und Ende der Lei- 
tung ist zum Vergleich in die Abb. 44 eingezeichnet. Die Zunahme 
der Dämpfung durch den erhöhten Leitungswiderstand infolge der 
Stromverdrängung ist nicht zu übersehen. 


Die Kurven wurden auf einem HP 67 mit dem folgenden Programm ge- 


rechnet 

Hauptprogranmn Ausgabe Unterpro- Unterprogramme 
Strom, gramm f. die Ränder 
Spannung BB i=0; i=4 

IBL A DSZ IASTX IBL O0 LBL c IBL a RCL 1 

9 STO 7 RCLA RC I POS RCL 2 - 

ST I - RI 2 RCL(i) RCL 2 RTN 

RCLE GSB 3 x A DSZ + 

1 GSB b LASTX PAUSE 3 RCLD 

+ IBL 4 R4 Rt} STO:(i) RCL O LBL b 

STO E RCL(i) + RCLB RCL(i) 2 RCL 9 

PAUSE 1SZ xOy : 8 x RCL 6 

2 RCL(i) POS —- STox(i) RCLC - 

B xQ 1S2Z Rt xDy x 2 

FRAC 704) ISZ -X- Rt - x 

x= xOy ISZ RTN STO+(i) - 0) 

GSB c DS2Z RTN LBL 3 DSZ LASTX GSB O0 

(6) SPACE LBL 2 RC I 3 xQy RTN 

ST I sto(i) STO(i) 2 STO:(i) RCLA 

GSB a RC I cLX B RCL(i) RCLC 

GSB O0 DSZ RCL A PAUSE STO+(i) - 

152 DSZ 1 R4 R4 RCL 1 

ISZ x£o0 + RCLB STO+(i) x 

GSB 1 GTo 4 STox(i) : R4 + 

GSB 2 GTOA Rh —- ISZ RCL A 

STO 3 LBL 1 DSZ DSZ STO+(i) 2 

- ISZ RCL(i) DSZ 4 + 

GSB 3 RCL(i 2 SPACE STO:(i) RCL C 

GSB 1 RCL(i x RCL(i) FO + 

ISZ RCL(i ISZ CHS DS2Z ; 

GSB 2 + STO-(i) xOI DS2 RCL 1 

DSZ 152 - 4 GTO c Rf 

STO 5 RCL(i) RCLA + RTN RCL 1 

- xSI 2 xQö Rt 

GSB 3 4 + RCL(i) STO 1 

GSB 1 - STO:(i) + + 

132 xG1l : —X- RCLC 

152 RCL(i) POS RTN x 

GSB 2 + DSZ + 

DS2 - REN Rt 


Das Zeittaktverhältnis ist der Programmschritt Nr. 9. Beim Vor- 
anschreiten von der Zeitzeile j zu j+1 müssen die Datenspeicher 
im Wanderwellengitter umgesetzt werden. Wegen der besonderen 
Aufteilung der Datenspeicher auf zwei Gittersysteme erfolgt 
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dieses in einem Unterprogramm: Schritt Nr. 42 bis 56. Das Pro- 
gramm enthält gelegentlich den Leerbefehl "SPACE", um einen 
Sprung nach DSZ auszugleichen. 


Speicherplan 
A B C D E I 
(34-1) Z, q E Zähler j Zähler i 
- 0,8755 50 1,6 10 Beginn: -1 


Die übrige Speichereinteilung entspricht der Abb. 53 


11 14 17 
12 15 18 
13 16 19 


0 2 4 6 8 
1 3 5 7 9 


10. Membranschwingung 


10.1 Problemstellung und Lösungsweg 


Eine am Rande eingespannte, dünne, elastische Haut (Membran) un- 


terliege einem einseitigen Druck, Sie wölbt sich 
einzelnen Membranpunkte erfahren eine Auslenkung 
die Membran plötzlich entlastet, so schwingt sie 
zeitliche Verlauf dieser Schwingung ist gesucht. 


aus, d.h. die 
ulx,y) . Wird 
zurück. Der 


Die Fragestellung führt auf die in der y- Richtung erweiterte 


Wellengleichung 
Ben. Du) 2.4250 
dx° dy° dt? 


mit Dichte £ und Spannung der Membran o: 


Beschreibung des Verfahrens 


Mit 
x=ih 
y=i'h 
= j-k 


(10.1) 


(10.2) 


werden die Differentialquotienten in der Gleichung (10.1) durch 


Zentralformeln ersetzt. Die Abmessungen der Gitterzelle h und k 
sollen dabei nicht unabhängig voneinander sein. Die Beziehung 


a (ie) 


(10.3) 
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führt auf die Rekursionsformel 


+ + 


) + 


= n-(u, 


%,i',jet 341,30,37% ea, ae ri 


+ 2-(1 - 2n)u. . (10.4) 


i,i' AP z ke Pe „j-1 
Die Rechnung ist für (1 - 2n)> O0 nmunmerisch stabil. Wählt man 
das Gleichheitszeichen, d.h. n 1/2 ,„ so vereinfacht sich die 
Rekursion: 


1 


u, ir,j+1 ” 2 ,ar,; bu w 


Mitar, Miitelss 


ringe 


Die Lösung ist eine unharmonische Schwingung. 


Programmablaufplan 


Die Membran habe einen rechteckigen Querschnitt a-b . Zu Beginn 
der Rechnung (Abb.55) müssen die Anfangswerte der beiden Zeit- 
schichten j = O0 und j = 1 eingegeben werden. Wie im Kapitel 6 
beschrieben, lassen sich die Werte durch "Auslenken" und "Los- 
lassen" finden. Um die Auslenkung u der einzelnen Membranpunkte 
für die Zeitschicht j = 0 zu erhalten, setzen wir die Membran 
unter einem einseitigen Druck p . Die Lösung der Poissonschen 
Differentialgleichung für diesen Fall nach Kapitel 11 ergibt 
die gesuchten u- Werte. Entlasten wir die Membran zur Zeit t# = 0, 
so ist die Rückschwing-Geschwindigkeit aller Punkte im ersten 
Moment noch Null. Dann ergeben sich, wie im Kapitel 6 beschrie- 
ben, die u- Werte der Zeitschicht j = 1 zu 


1 
Ur, = 2 Wapr,ir,o + Wiı,ir,o + %,ir41,0 + W,3r-1,0)° 


(10.6) 


Im weiteren Rechnungsablauf sind alle Ranäwerte Null zu setzen, 
weil die Membran an den Rändern eingespannt ist. Die u- Werte 
zweier aufeinanderfolgender j- Schichten werden zwischengespei- 
chert. Die Datenspeicher der Schicht j-1 dürfen dabei erst durch 
die Daten der Schicht j+1 überschrieben werden. 


-—4 speichern 


4; it,j+1 


Erhöhe 
um 1 


i' 


Erhöhe 
j um 1 


Abb. 55 Membranschwingung 
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10.2 Übergang auf Polarkoordinaten 


Ist die Membran kreisförmig, so überführen wir die Differential- 
gleichung nach Abschnitt 8.2 in Polarkoordinaten: 


2 2 
1 gu 2.09% u 
+. 5° = 5°; (10.7) 
r? dp? 


9 u BE RR 
Ir: ror 


Bei gleichmäßiger Belastung der Membran bleibt u unabhängig von 
y oder 


r=i.h 

t=j-k (10.8) 
2 

u 


zu der neuen Rekursion 


4 
4541 * n-[u4,1,; U Zi lauı,g” ] + (10,9) 


+ 2:(I-mu, ; u 031 


Diese ist numerisch stabil für (1-n)>0. Fürn= 1 erhält man 


5 RE 
Yo gti +t2 (4,5 - 41,5) = u ;_1° (10.10) 


Die Auslenkung u muß im Kreismittelpunkt (i = 0) endlich bleiben. 
Dazu ist 


%or1,j” Mo-1,j ” 
erforderlich, und wir erhalten im Kreismittelpunkt die Rekursion 


Uo,j+1 = ei: — Un ,j-1 . (10.11) 
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speichern 


speichern 


1,5 
Dur 
Y,j- -1 


Erhöhe 
jun 1 


Abb. 56 Membranschwingung, Übergang auf Polarkoordinaten 
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Der Programmablaufplan vereinfacht sich zur Abb. 56. Bei einem 
gleichmäßigen Druck auf eine Membran vom Außenradius a ergeben 
sich die u- Werte für die Zeitzeile j = O zu 


us=sc-(a-r) (10.12) 
und bei plötzlicher Entlastung für die Zeitzeile j = 1 


BE L. 
i,1=72 ar; +W1,0 + 21 (Urı1,o ” Be (10.13) 


und im Kreismittelpunkt 


Yo (10.14) 


10.3 Beispiel. Schwingung einer quadratischen Membran 


Die Schwingung einer quadratischen Membran, Kantenlänge a = 6 cm, 


b= 2,4103 s/cm „ soll berechnet werden. Wir unterteilen die 
Kante in sechs Teile. h = 1 cm ergibt nach (10.3) für n = 1/2 


k = 1,7-10”3 s . Eine quadratische Membran ist symmetrisch zum 

Achsenkreuz, so daß es genügt, die Rechnung für einen Quadranten 
auszuführen. Dieses erfordert neun Datenspeicher pro Zeitschicht 
(Abb.57). Dann müssen noch die Auslenkungen u der Zeitschichten 


y 


Abb. 57 Datenspeicher-Gitter bei einer quadratischen Membran 
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3j=-0 und j = I bestimmt werden. Hierbei sind wir in der “ahl 
des Auslenkungsdruckes frei. Mit dem Wert 


1 
f(x,y) = - = 0,39: 

und einer Genauigkeitsschranke von drei Stellen nach dem Komma 

ergibt das Gleichungssysten (11.4) folgende u- Werte für die 


Zeitschicht j = 0 in mm: 


- 0,597  - 0,545  - 0,370 
- 0,908  - 0,824  - 0,545 
- 1,003  - 0,908  - 0,597 . 


Die Entlastung nach Gleichung (10.6) liefert die u- Werte für 
die Zeitschicht j = I in mm: 


- 0,499 - 0,448 - 0,272 
=.0,812 - 0,727  - 0,448 
- 0,908 - 0,812 - 0,499 . 
u 
mm 
1,0 
0,5 
10 . 
J 
ke 0,5 
- 1,0 


Abb. 58 Schwingung des Membranmittelpunktes 
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Damit kann die Schwingung nach (10.5) gerechnet werden. Das Er- 
gebnis für den Membranmittelpunkt zeigt die Abb. 58. Als Schwin- 
gungszeit entnehmen wir Im = 12 Zeitschritte und erhalten damit 


die Frequenz nach (6,5) zu f=49 Hz. 


Gerechnet wurde auf einem HP 67. Die Ausgabe im Hauptprogramm 
beschränkte sich auf die Werte für den Membranmittelpunkt. Es 
ist aber nicht ohne Reiz, auch die Schwingungen der anderen Git- 
terpunkte mit in die Betrachtung einzubeziehen. 


Hauptprogramm Belastung EN IaBL ONE 
j=0 j = 

LBL A RCL 8 LBL B RCL 4 LBL C RCL 9 

RCLE 0 RCL 2 GSB 6 1 RCL 7 

1 RCL 4 RCL 2 STOo 7 1 (6) 

+ GSB 1 + RCL 9 ST I RCL 5 

STOE RBCL 9 RCL 4 + RCL 2 GSB O 

PAUSE RCL 7 RCL 4 [6) RCL 2 (6) 

1 16) GSB 6 RCL 5 RCL 4 RCL 8 

ST I RCL 5 STO 1 GSB 6 RCL 4 0 

RCL2 GsB 1 RCL 3 sto 8 GSB O0 RCL 6 

RCL 2 (6) + oO RCL 3 GSB 0 

RCL 4 RCL 8 RCL5 RCL 6 ROL 1 Pos 

RCL 4 (6) RCL 5 GSB 6 RCL 5 RTN 

GSB 1 RCL 6 GSB 6 RCL 9 RCL 5 

RCL 3 GSB 1 STO 2 ıOYy GSB O0 LBL O0 

RCL 1 POS RCL 6 STO 9 (6) + 

RCL 5 GTO A RCL 6 - RCL 2 + 

RCL 5 IBL 1 GSB 6 RND RCL 6 + 

GSB 1 + STO 3 x£o0 RCL 6 4 

[6) + RCL 5 GTO B GSB O0 : 

RCL 2 + RCL5 RTN RCL5 STO(i) 

RCL 6 2 + RCL 5 ISZ 

RCL 6 : RCL 7 LBL 6 RCL 7 RTN 

GSB 1 POS RCL 1 + RCL 1 

RCL 5 RCL(i) GSB 6 + GSB O0 

RCL 5 - STO 4 RCL A RCL 6 

RCL 7 STO(i) RCL 6 RCLD RCL 4 

RCL 1 GsB 2 + x RCL 8 

GSB 1 ISZ RCL 8 - RCL 2 

RCL 6 POS RCL 2 4 GSB 0 

RCL 4 RTN GSB 6 : (6) 

RCL 8 LBL 2 STO 5 RTN RCL 5 

RCL 2 RC I RCL 9 RCL 9 

GSB 1 1 RCL 3 RCL 3 

{6} xfYy GSB 6 GSB 0 

RCL 5 RTN STO 6 RCL 8 

RCL 9 PAUSE RCL 8 RCL 8 

RCL 3 RCL(i) RCL 8 (6) 

GSB 1 —- + RCL 4 

RCL 8 RTN 0 GSB O0 


Das Voranschreiten von j zu j+1 in den Datenspeichern erfolgt 
durch den Austausch der Primär- und Sekundärspeicher: POS . 
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Speicherplen 
A B [6 D E I 
frei frei h zähler j zöhler i,i' 
0,39 1 Beginn: 1 


Die übrige Speichereinteilung entspricht der Abb. 57. 


Zeitschicht j Zeitschicht j+1 
7 8 9 17 18 19 
4 5 6 14 15 16 
1 2 3 11 12 13 


Der Membranmittelpunkt ist der Datenspeicher 1 bzw. 11. 


11. Poissonsche- und Potentialgleichung 


11.1 Problemstellung und Lösungsweg 


Die Differentialgleichung der Wärmeleitung in einer dünnen Platte 
wurde im Kapitel 8 vorgestellt: 


2 2 
9°u u _ [) 
+ — = 2. , 


eg 


Es wird der Platte eine Temperaturverteilung u am Rande aufge- 
drückt. Nach einer gewissen Zeit hat sich im Innern eine statio- 
näre Temperaturverteilung eingestellt, die Ableitung nach der 
Zeit ist zu Null geworden: 


2 2 
au Tu 
+ —_.= [0] . (11.1) 
dx dy° 


Die Betrachtung auch anderer stationärer Felder mit einem Poten- 
tial u führt auf die Gleichung (11.1), die den Namen Potential- 
gleichung erhalten hat. Sie gilt aber nur, wenn im Innern der 
Platte keine Wärme erzeugt oder abgeführt wird. Gibt es im Innert 
einen Wärmefluß, so ist die rechte Seite der Gleichung nicht mehı 
Null: 


2 2 
gu ou 
+ —— = f(x,y) . (11.2) 
Ix° dy° : 


Die Funktion f(x,y) drückt dann den stationären Wärmefluß in den 
Plattenpunkten aus. Die um f(x,y) erweiterte Potentialgleichung 
heißt Poissonsche Differentialgleichung. Beide Gleichungen spie- 
len in vielen Anwendungsgebieten eine fundamentale Rolle. So 
führt die Druckbelastung einer Membran ebenfalls auf die Glei- 


11.1 Problemstellung und Lösungswer = 


chung (11.2). Mit dem Druck p, der Membrandicke d, der Spannurz : 
ist 


f(x,y) = - Ir 


und u(x,y) die Auslenkung der Membranpunkte. 
Für eine kreisrunde Umrandung schreibt sich (11.2) in Polarkoor- 
dinaten 


au 1 du BON = 
HE 5 Be 2 or (11.3) 


o 
ie} 


Ist u unabhängig von p „ so geht die partielle in eine gewöhn- 
liche Differentialgleichung zweiter Ordnung mit Randwerten über. 


Beschreibung des Verfahrens 


Es genügt, die Poissonsche Differentialgleichung zu behandeln, 
weil die Potentialgleichung in ihr als Sonderfall mit f(x,y) = 0 
enthalten ist. Die Differentialquotienten werden durch Zentral- 
formeln ersetzt. Die Abmessungen der Gitterzelle h und k sollen 
von einander unabhängig wählbar sein. Mit 


x=i-h 


y=jk 


erhält man eine Gleichung für u, in Gitterpunkt i,j: 


1,) 


2 
+ 44,5) +h (we + u, 5-1) + 


u 2 Be 
2-(h’+ k Ju; 5 =k (w,j 


- neretla,y) - (11.4) 


Die Abb. 59 soll diese Gleichung veranschaulichen. Der u- Wert 
der linken Seite ist durch einen leeren Kreis, die u- Werte der 
rechten Seite sind durch volle Kreise gekennzeichnet. 

Schreibt man die Gleichungen für alle Gitterpunkte i,j unterein- 


ander, indem man jeweils einen anderen u. j” Wert auf der linken 


1, 
Seite hat, so ergeben sich n Gleichungen mit n Unbekannten. Die 
Anzahl der Gitterpunkte im Imnern des Gitters ist gleich n. Die 
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Abb. 59 Zur Gleichung (11.4) 


Funktionswerte u in den Rand-Gitterpunkten und die Funktion 
f(x,y) sind nach Voraussetzung gegeben. 

Ein solches umfangreiches Gleichungssystem läßt sich auf einem 
Taschenrechner nur durch ein Iterationsverfahren lösen. Wegen 
seines einfachen Aufbaues bietet sich das Einzelschrittverfahren 
nach Gauß-Seidel an. 

Auf den linken Seiten des Gleichungssystems steht jeweils eine 
andere Unbekannte u. Der Startwert für alle Unbekannten ist Null. 
Mit den vorhandenen Zahlenwerten in den Rand-Gitterpunkten werden 
für die einzelner Unbekannten Näherungswerte, Gleichung hinter 
Gleichung, ermittelt. Hierbei sind stets die neusten Näherungs- 
werte in die nachfolgenden Gleichungen einzusetzen. Ist man auf 
diesem Wege bei der letzten Gleichung des Systems angekommen, so 
beginnt die Rechnung mit der ersten wieder von vorne. Die Rech- 
nung ist beendet, wenn sich, gemessen an einer Genauigkeitsschran 
ke, die Zahlenwerte der Unbekannter nicht mehr ändern. 

Eine schnelle Konvergenz dieses einfachen Iterationsverfahrens 
ist nicht bei jedem Gleichungssystem gewährleistet. Erforderlich 
sind wenige Unbekannte pro Gleichung und ein Überwiegen des Koef- 
fizienten auf der linken Seite gegenüber denen auf der rechten. 
Dieses trifft bei dem vorliegenden Gleichungssystem (11.4) zu, 
wenn nicht zu viele Gitterpunkte weit ab vom Rande liegen. Als 
Maß zur Beendigung der Rechnung kann man den u- Wert der letzten 
Gleichung wählen. Erreicht der Unterschied aus zwei aufeinander- 
folgender Iterationen 


Im m,,;) (11.5) 


m+1 


A= u;, m 


J i1,J 
eine vorgegebene Genauigkeitsschranke, so wird die Rechnung be- 
endet. 
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h, k 
7] Randwerte 


Mi-1,5 9 Yirt,z 


u ,j-1 9 %,je1 


Entfällt bei 
7 Potential- 
gleichung 


speichern 


60 Poissonsche- und 
Potentialgleichung 
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Programmablaufplan 


Das Gitter habe die Abmessung a-b . Zu Beginn der Rechnung 
(Abb.60) werden die Werte u an den vier Rändern sowie h und k 
eingegeben. Bei der Poissonschen Differentialgleichung kann 


die Funktion f, j in den einzelnen Gitterpunkten in einem Un- 
’ 


terprogramm gerechnet werden. Bei der Fotentialgleichung ent- 
fällt dieses. Während jedes Iterationsdurchlaufes überschreiben 
die neuen Näherungswerte sofort die alten in den entsprechenden 
Datenspeichern, mit Ausnahme des letzten Gitterpunktes. Hier wird 
zunächst der Unterschied nach (11.5) festgestellt und mit der Ge- 
nauigkeitsschranke verglichen. 


11.2 Beispiel. Temperaturverteilung in einem Wohnraum im Winter 


Ein Beispiel für die Poissonsche Differentialgleichung war die 
Auswölbung einer Membran innerhalb des Beispieles 10.3. Als Bei- 
spiel für die Potentialgleichung soll der Einfluß einer dünnen 
Balkontür auf die Temperaturverteilung in einem geheizten Wohn- 
raum ermittelt werden. Bekannt sind die Randtemperaturen, der 
niedrigste Wert an der Balkontür, der höchste an der Heizung. 
Auch weitere Gegebenheiten, wie Schornstein und Küche, werden in 
die Überlegungen einbezogen. 
Das Zimmer wird mit einem Gitter aus quadratischen Zellen über- 
zogen. Ohne die Ränder hat es 18 Punkte (Abb.61). Für h = k und 
f(x,y) = 0 vereinfacht sich (11.4) 

4-u +Uu 


EUR Banle 5 u Pr En Zn 5 Pr zu zn 92 22 Enz 9 u Be 
Gerechnet wird das System aus 18 Gleichungen nach dem beschrie- 
benen Iterationsverfahren mit einer Genauigkeitsschranke von 
zwei Stellen nach dem Komma. Das Ergebnis ist in Abb. 61 festge- 
halten. Die kalte Balkontür macht sich im Zimmer nur wenig be- 
merkbar, Dafür ist der Wärmestau in der Zimmerecke an der Hei- 
zung nicht zu übersehen, Der Einbau eines Thermostaten an dieser 
Stelle ist nicht empfehlenswert. 
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26 20 20 20 20 


Schornstein 


32 

Küche Heizung 
32 
19 

Flur Balkon 


9 


20 20 20 20 


Abb. 61 Datenspeicher-Gitter und Temperaturverteilung 
des Beispiels 11.2 


Gerechnet wurde die Temperaturverteilung auf einem HP 67 mit dem 
folgenden Iterationsprogram für 18 Gleichungen mit 18 Unbekann- 
ten. Die Genauigkeitsschranke entspricht dem Anzeigeformat. 


144 11. Poissonsche- und Potentialgleichung 


LIBL A GSB 1 GSB 1 sto 0 RC I RCLD 
RCL 2 STO 4 STo 7 + + 
RCLB RcL 2 RCLE PO 
+ POs RCLC + RCL 5 RCL A 
RCLA RCL 2 + RCL 6 GSB 1 RCL 9 
RCL 8 POS RCL 9 POS STO 4 GSB 1 
GSB 1 + RCL 1 RCL 6 PQS 
STO 1 RCLA GSB 1 FO RCL 6 RCL 8 
RCL 6 sTo 8 GSB 1 + xy 
RCL 3 GSB 1 sTto 1 RCLA sto 8 
+ STo 5 POS RCL 2 PQSs 
RcL 7 RCL 8 RCL 3 GSB 1 Be 
RCLA RCL 7 + + STO 5 RND 
GSB 1 + RCLD RCL 5 x#£0 
STo 2 Pos RCLO POs RCL 7 GTO A 
RCL 1 GSB 1 RCL 5 + RTN 
RCL 4 püs POS PS RCLA 
+ RCL 3 sto 9 GSB 1 RCL 1 
RCL 6 GSsB 1 sTto 2 GSB 1 LBL 1 
RCLA STO 6 POS STO 6 + 
GSB 1 RCL 1 RCLE + 
sto 3 RCL 8 + + RCL 8 4 
+ RCL 7 RCL 4 + : 
RCLA Pos POS RCLA RCLA RTN 
+ RCL Oo RCL 7 GSB 1 RCL O 
RCL 5 POS PGs STO 3 GSB 1 
RCLA RCL 2 GSB 1 STo 7 
Speicherplan 
A B C D E I 
Randwerte 20 9 19 32 23 26 


Die übrige Speichereinteilung entspricht der Abb. 61. 


14 15 16 17 18 
13 12 11 10 

5 6 7 

4 3 2 1 
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